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Chapter 0 

outline of this lecture 


This document is the lecture script of a one-semester course taught at the University 
of Basel in the Fall semesters of 2012 and 2013 and in the Spring semester of 2015. 
It is aimed at advanced students of physics who are familiar with the concepts and 
notations of quantum mechanics. 


0.1 introduction 

Quantum mechanics lectures can often be separated into two classes. In the first 
class you get to know Schrodinger’s equation and find the form and dynamics of 
simple physical systems (square well, harmonic oscillator, hydrogen atom); most 
calculations are analytic and inspired by calculations originally done in the 1920s 
and 1930s. In the second class you learn about large systems such as molecular 
structures, crystalline solids, or lattice models; these calculations are usually so com- 
plicated that it is difficult for the student to understand them in all detail. 

This lecture tries to bridge the gap between simple analytic calculations and 
brute-force large-scale computations. We will revisit most of the problems encoun- 
tered in introductory quantum mechanics, focusing on computer implementations 
for finding analytical as well as numerical solutions and their visualization. We will 
be approaching topics such as the following: 

• You have calculated the energy eigenstates of single particles in simple poten- 
tials. How can such calculations be generalized to non-trivial potentials? 

• How can we calculate the behavior of interacting particles? 

• How can we describe the internal spin structure of particles? How does this 
internal structure couple to the particles’ motion? 

• You have heard that quantum mechanics describes our everyday world just 
as well as classical mechanics does, but you may never have seen an example 
where this is calculated in detail and where the transition from the classical 
behavior to the quantum-mechanical behavior is evident. 

Most of these calculations are too complicated to be done by hand. Even rela- 
tively simple problems, such as two interacting particles in a one-dimensional trap, 
do not have analytic solutions and require the use of computers for their solution 
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and visualization. More complex problems scale exponentially with the number of 
degrees of freedom, and make the use of large computer simulations unavoidable. 


0.2 what this lecture is not about 

This course is not about quantum computing. Quantum computing refers to the 
proposed use of quantum-mechanical entanglement of physical systems for speed- 
ing up certain calculations, such as factoring large numbers. 

This course is not about large-scale quantum calculations such as solid-state 
physics or quantum chemistry. It will, however, provide you with the tools for un- 
derstanding such large-scale calculations better. 


0.3 Why Mathematica? 

The course will be taught in the Wolfram language of Mathematica (version 10) . No 
prior knowledge of Mathematica is necessary, and Mathematica licenses will be pro- 
vided. Alternatives to Mathematica, such as Matlab or Maple, may be used by the 
students, but only limited help will be available from the instructor. 

There are many reasons for choosing Mathematica over other computer-algebra 
systems (CAS): 

• Mathematica is a very high-level programming environment, which allows the 
user to focus on what he wants to do instead of how it is done. A very large 
number of algorithms for analytic and numerical calculations is included in 
the Mathematica kernel and its libraries. 

• Mathematica seamlessly mixes analytic and numerical facilities. For many 
calculations it allows you to push analytic evaluations as far as possible, and 
then continue with numerical evaluations by making only trivial changes. 

• Mathematica supports a wide range of programming paradigms, which means 
that you can keep programming in your favorite style. See subsection 1.6.5 for 
a concrete example. 

• The instructor is more familiar with Mathematica than any other CAS. 

0.4 outline of discussed topics 

Chapter 1 gives an introduction to Mathematica and the Wolfram language, with a 
focus on techniques that will be useful for this lecture. 

Chapter 2 makes the connection between quantum mechanics and the vector/matrix 
calculus of Mathematica. 

Chapter 3 discusses systems with finite-dimensional Hilbert spaces, focusing on 
spin systems. 

Chapter 4 discusses the quantum mechanics of particles moving in one- and several- 
dimensional space. 

Chapter 5 connects the topics of chapters 3 and 4. 


0 . 4 . OUTLINE OF DISCUSSED TOPICS 
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Chapter 6 presents a brief introduction to path integrals and Monte Carlo integra- 
tion. 
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Chapter 1 

Wolfram language overview 


If you have little or no experience with Mathematica and the Wolfram language, you 
may start here: 

http : / /www . wolfram. com/support/learn/higher- education .html 
Further useful links: 

• http : //www . wolfram . com/language/ 

• http : //reference . wolfram. com/mathematica/guide/LanguageOverview . 
html 

• http : //reference . wolfram. com/mathematica/guide/Mathematica. html 


1.1 introduction 

Mathematica is an interactive system for mathematical calculations. The Mathe- 
matica system is composed of two main components: the front end, where you write 
the input in the Wolfram language, give execution commands, and see the output, 
and the kernel, which does the actual calculations. 


front end kernel 



This distinction is important to remember because the kernel remembers all the 
operations in the order they are sent to it, and this order may have nothing to do 
with the order in which these commands are displayed in the front end. 

When you start Mathematica you see an empty “notebook” in which you can 
write commands. These commands are written in a mixture of text and mathemat- 
ical symbols and structures, and it takes a bit of practice to master all the special 
input commands. In the beginning you can write all your input in pure text mode, if 
you prefer. Let’s try an example: add the numbers 2 + 3 by giving the input 


11 


12 


CHAPTER 1 . WOLFRAM LANGUAGE OVERVIEW 


In[l] := 2+3 


and, with the cursor anywhere within the "cell” containing this text (look on the right 
edge of the notebook to see cell limits and groupings) you press “shift- enter”. This 
sends the contents of this cell to the kernel, which executes it and returns a result 
that is displayed in the next cell: 


Qut[i]= 5 


If there are many input cells in a notebook, they only get executed in order if you 
select “Evaluate Notebook” from the “Evaluation” menu; otherwise you can execute 
the input cells in any order you wish by simply setting the cursor within one cell and 
pressing “shift- enter”. 

1.1.1 exercises 


Do the following calculations in Mathematica, and try to understand their structure: 
Q 1 . 1 Calculate the numerical value of ( (3) with 



Q1.3 Calculate sin(x)e ^dxwith 


i in [4] : = Integrate [Sin [x] Exp[-x], {x, 0, Infinity}] 


Q1.4 Calculate the first 1000 digits of n with 


i In [5] := N [Pi , 1000] 


Q1.5 Calculate the Clebsch-Gordan coefficient <1000, 100; 2000, — 120| 1 100, -20): 


i in [6] := ClebschGordan [{1000, 100}, {2000, -120}, {1100, -20}] 


Q1.6 Calculate the limit lim ,_o — p- with 


i in [7] := Limit [Sin [x] /x, x -> 0] 


Q1.7 Make a plot of the above function with 


i in [8] := Plot [Sin [x] /x , {x , -20, 20}, PlotRange -> All] 
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Q1.8 Draw a Mandelbrot set with 

i in [9] := F[c_, imax_] := Abs [NestWhile [#~2 + c k, 0., 
z Abs [#] <= 2 &, 1, imax]] <= 2 

3 in [io] : = With[{n = 100, imax = 1000}, 

4 Graphics [Raster [Table [Boole [ !F [x + I y, imax]], 

5 {y, -2, 2, 1/n}, {x, -2, 2, 1/n}]]]] 


Q1.9 Do the same with a built-in function call: 


i in[n] := MandelbrotSetPlot [] 


1 .2 variables and assignments 


http : / /reference . wolfram . com/mathematica/howto/WorkWithVariablesAndFunctions . html 

Variables in Mathematica can be letters or words with uppercase or lowercase 
letters, including Greek symbols. Assigning a value to a variable is done with the = 
symbol, 

In [12] := a = 5 
Out [12]= 5 


If you wish to suppress the output, then you must end the command with a semi- 
colon: 

In [13] := a = 5; 


The variable name can then be used anywhere in an expression: 

In [14]:= a + 2 
Out [14]= 7 


1.2.1 immediate vs. delayed assignments 

http : / / reference . wolfram . com/mathematica/tutorial/ImmediateAndDelayedDef init ions . html 

Consider the two commands 

in [ 15 ] := a = RandomRealf] 

Out [ 15 ]= 0.38953 

in[i6] := b := RandomRealf] 


(your random number will be different). 

The first statement a= . . . is an immediate assignment, which means that its 
right-hand side is evaluated when you press shift-enter, produces a specific random 
value, and is assigned to the variable a (and printed out). From now on, every time 
you use the variable a, the exact same number will be substituted. In this sense, 
the variable a contains the number 0.38953 and has no memory of where it got this 
number from. You can check the definition of a with ?a: 
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In [17] := ?a 

Global ‘a 
a = 0.38953 


The definition b : = . . . is a delayed assignment, which means that when you 
press shift-enter the right-hand side is not evaluated but merely stored as a defi- 
nition of b. From now on, every time you use the variable b, its right-hand-side 
definition will be substituted and executed, resulting in a new random number each 
time. You can check the definition of b with 

In [18]:= ?b 

Global ‘b 

b : = RandomReal [] 


Let’s compare the repeated performance of a and b: 


In [19] := fa, b> 

Out [ 19 ] = -[0.38953, 0.76226} 

In [20]:= {a, b} 

Out [ 20 ] = -[0.38953, 0.982921} 
In [21] := {a, b} 

Out [ 21 ] = {0.38953, 0.516703} 

In [22]:= {a, b} 

Out [ 22 ] = {0.38953, 0.0865169} 


1.2.2 exercises 

Q1.10 Explain the difference between 



In particular, distinguish the cases where u and v are already defined before x 
and y are defined, where they are defined only afterwards, and where they are 
debited before but change values after the debnition of x and y. 

1 .3 four kinds of bracketing 

http : / / reference . wolfram . com/language/tutor ial/TheFourKindsOf Bracket inglnTheWolf ramLanguage . html 

There are four types of brackets in Mathematica: 

• parentheses for grouping, for example in mathematical expressions: 

i In [25] := 2*(3-7 ) 
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• square brackets for function calls: 
i In [ 26 ] := Sin [0.2] 


• double square brackets for indexing within lists: (see section 1.7) 

i In [28] := V [ [2] ] 


1 .4 prefix and postfix 

There are several ways of evaluating a function call in Mathematica, and we will see 
most of them in this lecture. As examples of function calls with a single argument, 
the main ways in which sin(0.2) and v2 + 3 can be calculated are 

standard notation (infinite precedence): 

1 
2 

3 

4 


In [29] := Sin [0.2] 
out [ 29 ]= 0.198669 
in [ 30 ] := Sqrt [2+3] 
out [30]= Sqrt [5] 


curly braces for lists: 


In [27] := V = {a, b, c} 


prefix notation with @ (quite high precedence, higher than multiplication): 


2 

3 

4 


Notice how the high precedence of the @ operator effectively evaluates ( Sqrt @2 ) +3, 
not Sqrt® (2+3). 

postfix notation with // (quite low precedence, lower than addition): 


2 

3 

4 


In [33] := 0.2 // Sin 
Out [ 33 ] = 0.198669 
In [34] := 2+3 // Sqrt 
out [34]= Sqrt [5] 


In[31] := Sin @0.2 
Out [ 31 ] = 0.198669 
in [ 32 ] := Sqrt @ 2+3 
out [32]= 3+Sqrt [2] 


Notice how the low precedence of the // operator effectively evaluates (2+3) //N, 
not 2+(3//N). 

Postfix notation is often used to transform the output of a calculation: 

• Adding //N to the end of a command will convert the result to decimal 
representation, if possible. 
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• Adding //Matr ixForm to the end of a matrix calculation will display the 
matrix in a tabular form. 

• Adding //Timing to the end of a calculation will display the result to- 
gether with the amount of time it took to execute. 

If you are not sure which form is appropriate, for example if you don’t know the 
precedence of the involved operations, then you should use the standard notation 
or place parentheses where needed. 

1.4.1 exercises 

Ql.ll Calculate the decimal value of Euler’s constant e (E) using standard, prefix, 
and postfix notation. 


1 . 5 programming constructs 

When you program in Mathematica you can choose between a number of differ- 
ent programming paradigms, and you can mix these as you like. Depending on the 
chosen style, your program may run much faster or much slower. 

1.5.1 procedural programming 

http : / / reference . wolfram . com/ mathematica/ guide/ProceduralProgramming . html 

A subset of Mathematica behaves very similarly to C, Python, Java, or other pro- 
cedural programming languages. Be very careful to distinguish semi-colons, which 
separate commands within a single block of code, from commas, which separate 
different code blocks! 

Looping constructs behave like in common programming languages: 



Conditioned execution: notice that the If statement has a return value, similar to 
the “?” statement of C and Java. 


2 

3 


In [38] : = If [5! > 100, 

Print ["larger"] , 

Print ["smaller or equal"]] 


1.5. PROGRAMMING CONSTRUCTS 
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2 


Modularity: code can use local variables within a module'. 


2 

3 

4 

5 


In [40]:= Module [{i}, 
i = 1; 

While [i > 1/192, i = i/2]; 
i] 

out [ 40 ]= 1/256 


In [39] := a = If [5! > 100, 1, -1] 

Out [39]= 1 


After the execution of this code, the variable i is still undefined in the global 
context. 

1.5.2 exercises 

Q1.12 Write a program which sums all integers from 123 to 9968. Use only local vari- 
ables. 

Q1.13 Write a program which sums consecutive integers, starting from 123, until the 
sum is larger than 10000. Return the largest integer in this sum. Use only local 
variables. 

1.5.3 functional programming 

http : / /reference . wolfram . com/mathematica/ guide/Funct ionalProgramming . html 

Functional programming is a very powerful programming technique which can 
give large speedups in computation because it can often be parallelized over many 
computers or CPUs. In our context, we often use lists (vectors or matrices, see sec- 
tion 1.7) and want to apply functions to each one of their elements. 

The most common functional programming constructs are 

Anonymous functions: 1 you can quickly define a function with parameters #1=#, 
#2, #3, etc., terminated with the k symbol: 



Functions and anonymous functions, for example #"2&, are first-class ob- 
jects * 2 just like numbers, matrices, etc. You can assign them to variables, as 
above; you can also use them directly as arguments to other functions, as be- 
low. 

The symbol ## stands for the sequence of all parameters of a function: 


x see http : //en. Wikipedia. org/wiki/ Anonymous .functions. 

2 see http : //en. Wikipedia. org/wiki/Fir st- class.citizen. 
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1 In [45] : = f = { 1 , 2 , 3 , ## , 4 , 5 , 6} &; 

2 In [46] : = f [7 , a, c] 

3 out [46]= {1, 2, 3,7, a, c, 4, 5, 65- 


Map /<§: apply a function to each element of a list. 


In [47] : = 

a = 

{1, 2, 3, 

4, 

5, 6, 

7, 

8>; 

In [48] : = 

Map 

1 — 1 
aJ 

CM 

< 





Out [48] = 

{1, 

4, 9, 16, 

25, 

36, 

49, 

64} 

In [49] : = 

#~2 

& /@ a 





Out [49] = 

{1, 

4, 9, 16, 

25, 

36, 

49, 

64} 


Notice how we have used the anonymous function #~2 k here without ever 
giving it a name. 

Apply @0: apply a function to an entire list and generate a single result. For example, 
applying Plus to a list will calculate the sum of the list elements; applying 
Times will calculate their product. This operation is also known as reduce:’ 


1 

in [ 50 ]:= a = {1, 2, 3, 4, 5, 6, 7, 8}; 

2 

In [51]:= Apply [Plus, a] 

3 

Out [5i]= 36 

4 

in [52]:= Plus @0 a 

5 

Out [52]= 36 

6 

in [53] : = Apply [Times, a] 

7 

out [53]= 40320 

8 

in [54] : = Times @@ a 

9 

Out [54]= 40320 


1.5.4 exercises 

Q1.14 Write an anonymous function with three arguments, which returns the prod- 
uct of these arguments. 

Q1.15 Given a list 

1 in [55] : = a = {0.1, 0.9, 2.25, -1.9}; 

calculate x *— * sin(x 2 ) for each element of a using the Map operation. 

Q1.16 Calculate the sum of all the results of Q1.15. 


1.6 function definitions 


http: //ref erence .wolf ram. com/mathematica/tutorial/Def iningFunctions .html 

Functions are assignments (see section 1.2) with parameters. As for parameter- 
free assignments we distinguish between immediate and delayed function defini- 
tions. 

3 See http : / /en . Wikipedia . org/wiki/MapReduce. 
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1.6.1 immediate function definitions 

We start with immediate definitions: a function f{x ) = sin(x)/x is defined with 
In [56] := f [x_] = Sin [x] /x ; 


Notice the underscore _ symbol after the variable name x: this underscore indicates 
a pattern (denoted by _) named x, not the symbol x itself. Whenever this function 
f is called with any parameter value, this parameter value is inserted wherever x 
appears on the right-hand side, as is expected for a function definition. You can find 
out how f is defined with the ? operator: 


In [57] := ?f 

Global ‘f 

f [x_] = Sin[x]/x 


and you can ask for a function evaluation with 


In [58] := f [0 . 3] 

Out [58]= 0.985067 
In [59] := f [0] 

Power: : inf y : Infinite expression 1/0 encountered. 

Inf inity : : indet : Indeterminate expression 0 Complexlnf inity encountered. 

out [59]= Indeterminate 


Apparently the function cannot be evaluated for x - 0. We can fix this by defining a 
special function value: 

In [60] := f [0] = 1 J 


Notice that there is no underscore on the left-hand side, so there is no pattern defi- 
nition. The full definition of f is now 


In[61] := ?f 

Global ‘f 
f CO] =1 

f [x_] = Sin[x]/x 


If the function f is called, then these definitions are checked in order of appearance 
in this list. For example, if we ask for f [0] , then the first entry matches and the 
value 1 is returned. If we ask for f [0 . 3] , then the first entry does not match (since 
0 and 0.3 are not strictly equal), but the second entry matches since anything can 
be plugged into the pattern named x. The result is sin(0.3)/0.3 = 0.985067, which is 
what we expected. 

1.6.2 delayed function definitions 

Just like with delayed assignments (subsection 1.2.1), we can define delayed func- 
tion calls. For comparison, we define the two functions 
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in [ 62 ] := gl [x_] = x + RandomReal [] 

out [ 62 ]= 0.949868 + x 

in [63] : = g2 [x_] := x + RandomReal [] 


Check their effective definitions with ?gl and ?g2, and notice that the definition of 
gl was executed immediately when you pressed shift-enter and its result assigned 
to the function gl (with a specific value for the random number, as printed out), 
whereas the definition of g2 was left unevaluated and is executed each time anew 
when you use the function g2: 

In [64]:= {gl [2] , g2 [2] > 

Out [64] = {2.94987, 2.33811} 

In [65]:= {gl [2] , g2 [2] } 

Out [65] = {2.94987, 2.96273} 

In [66]:= {gl [2] , g2 [2] } 

Out [66] = {2.94987, 2.18215} 


1.6.3 functions that remember their results 

http : //ref erence .wolf ram. com/mathematica/tutorial/Funct ionsThatRememberValuesTheyHaveFound.html 

When we define a function that takes a long time to evaluate, we may wish to 
store its output values such that if the function is called with identical parameter 
values again, then we do not need to re-evaluate the function but can simply return 
the already calculated result. We can make use of the interplay between patterns 
and values, and between immediate and delayed assignments, to construct such a 
function which remembers its values from previous function calls. 

See if you can understand the following definition. 



and ask for ?F again. You see that the specific immediate definition of F [2] =128 
was added to the list of definitions, with the evaluated result 128 (which may have 
taken a long time to calculate in a more complicated function). The next time you 
call F [2] , the specific definition of F [2] will be found earlier in the definitions list 
than the general definition F [x_] and therefore the precomputed value of F [2] will 
be returned. 

When you re-define the function F after making modifications to it, you must 
clear the associated remembered values in order for them to be re-computed at the 
next occasion. It is a good practice to prefix every definition of a self-remembering 
function with a Clear command: 

in [69] := Clear [F] ; 

In [70] : = F[x_] := F [x] = x~9 


1.6. FUNCTION DEFINITIONS 
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1.6.4 functions with conditions on their arguments 

http : / / reference . wolfram . com/mathemat ica/ guide/Patterns . html 

Mathematica contains a powerful pattern language that we can use to define 
functions which only accept certain arguments. For function definitions we will use 
three main types of patterns: 

Anything-goes: A function defined as 


i In [71] := f [x_] : = X~2 


can be called with any sort of arguments, since the pattern x_ can match any- 
thing: 


1 

In [72]:= f [4] 

2 

Out [72]= 16 

3 

In [73] := f [2.3-0. 11] 

4 

Out [73] = 5.28-0.461 

5 

In [74] := f [{1,2, 3, 4}] 

6 

Out [74] = {1,4,9,16} 

7 

In [75] := f [y“2] 

8 

Out [75]= y~4 


Type-restricted: A pattern like x_Integer will match only arguments of integer 
type. If the function is called with a non-matching argument, then the func- 
tion is not executed: 


1 

in [76]:= g[x_Integer] := x-3 

2 

g[x_Real] := x+3 

3 

In [77] := g [7] 

4 

Out [77]= 4 

5 

In [78]:= g [7 . 1] 

6 

Out [78] = 10.1 

7 

In [79]:= g[2+3I] 

8 

Out [79]= g [2+31] 


Conditional: Complicated conditions can be specified with the / ; operator: 


In [80] : = 

h[x_/ 

x<=3] 

:= x~2 


h[x_/ 

1 1 

co 

A 

X 

= X- 1 1 

In[81] : = 

h [2] 



Out [81] = 

4 



In [82] : = 

h[5] 



Out [82] = 

-6 




Conditions involving a single function call returning a Boolean value, for ex- 
ample x_/ ; Pr imeQ [x] , can be abbreviated with x_?PrimeQ. 
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1.6.5 fifteen ways to define the factorial function 

The Wolfram demo http: //www . wolf ram. com/training/videos/EDU002/ lists 
fifteen ways how to define the factorial function. Try to understand as many of these 
definitions as possible. What this means in practice is that for most problems you 
can pick the programming paradigm which suits your way of thinking best, instead 
of being forced into one way or another. The different paradigms have different ad- 
vantages and disadvantages, which may become clearer to you as you become more 
familiar with them. 

You must call Clear [f ] between different definitions! 

1. Define the function f to be an alias ofthe built-in functionFactorial: calling 
f [5] is now strictly the same thing as calling Factorial [5] , which in turn is 
the same thing as calling 5 ! . 


i in [83] : = f = Factorial; 


2. A call to f is forwarded to the function “ ! calling f [5] triggers the evaluation 
of 5!. 

i in [84] := f [n_] := n! 


3. Use the mathematical definition n\ — T(n + 1): 
i in [85] : = f [n_] := Gamma [n+1] 


4. Use the mathematical definition n\ — n " =1 i: 
i in [86] : = f [n_] := Product [i , {i,n}-] 


5. Rule-based recursion, using Mathematica’s built-in pattern-matching capa- 
bilities: calling f [5] leads to a call of f [4] , which leads to a call of f [3] , and 
so on until f [1] immediately returns the result 1, after which the program 
unrolls the recursion stack and does the necessary multiplications: 


2 


6. The same recursion but without rules (no pattern-matching): 
i in [88]:= f [n_] := If [n == 1, 1, n f [n-1]] 


7. Define the same recursion defined through functional programming: f is a 
function whose name is #0 and whose first (and only) argument is #1. The 
end of the function definition is marked with &. 

i In [89]:= f = If [# 1 == 1, 1, #1 #0 C# 1 - 1 ] ] & J 


In [87] := f [1] = 1; 

f [n_] : = n f [n-1] 


1.6. FUNCTION DEFINITIONS 
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i. procedural programming with a Do loop: 


in [90] : = f [n_] := Module [{t = 1}, 
Do [t = t i, {i, n>] ; 
t] 

procedural programming with a For loop: this is how you would compute fac- 
torials in procedural programming languages like C. It is a very precise step- 
by-step prescription of how exactly the computer is supposed to do the calcu- 
lation. 

in [ 91 ] := f [n_] := Module [{t = 1, i}, 
For[i = 1 , i <= n, i++, 
t *= i] ; 
t] 

Make a list of the numbers 1 . . . n (with Range [n] ) and then multiply them to- 
gether at once, by applying the function Times to this list. This is the most 
elegant way of multiplying all these numbers together, because both the gen- 
eration of the list of integers and their multiplication are done with internally 
optimized methods. The programmer merely specifies what he would like the 
computer to do, and not how it is to be done. 

in [92] := f [n_] := Times @@ Range [n] 

Make a list of the numbers 1 . . . n and then multiply them together one after 
the other. 

in [93] := f [n_] := Fold [Times, 1, Range [n] ] 

Functional programming: make a list of functions {t<-~ t, t >-► 2 1 , t <->■ 3 1, . . . , t >-► 
n t}, and then, starting with the number 1, apply each of these functions once. 

in [94] := f [n_] := Fold [#2 [#1] & , 1, 

Array [Function [t , #1 t]&, n] ] 

Construct a list whose length we know to be n\: 

in [95] := f [n_] := Length [Permutat ions [Range [n] ] ] 

Use repeated pattern-based replacement (/ / . ) to find the factorial: start with 
the object {1, n) and apply the given rule until the result no longer changes 
because the pattern no longer matches. 

in [ 96 ]:= f [n_] := First [{l,n} // . {a_ ,b_/ ;b> 0 } :> {b a,b-l}] 


15. Build a string whose length is n\: 
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1 in [ 97 ] : = f [n_] := Str ingLength [Fold [ 

2 StringJoin [Table [#1 , {#2}]]&, 

3 "A", Range [n]]] 


1.6.6 exercises 

Q1.17 In which ones of the definitions of subsection 1.6.5 can you replace a de- 
layed assignment ( : =) with an immediate assignment (=) or vice-versa? What 
changes if you do this replacement? 

Q1.18 Can you use the trick of subsection 1.6.3 for any of the definitions of subsec- 
tion 1.6.5? 

Q1.19 Write two very different programs that calculate the first hundred Fibonacci 
numbers { 1 , 1 , 2 , 3, 5, 8 , . . .}, where each number is the sum of the two preceding 
ones. 

1 .7 vectors and matrices 

In this lecture we will use vectors and matrices to represent quantum states and 
operators, respectively. 

1.7.1 vectors 

http : / / reference . wolfram . com/mathemat ica/ tutorial/VectorOperat ions . html 

In Mathematica, vectors are represented as lists of objects, for example lists of 
real or complex numbers: 

In [98]:= V = {1,2, 3, 2, 1,7+1}-; 
in [ 99 ] := Length [v] 

Out [99] = 6 


You can access any element by its index, using double brackets, with the first ele- 
ment having index 1 (as in Fortran or Matlab), not 0 (as in C, Java, or Python): 

In[100] : = V [ [4] ] 
out [ioo] = 2 


Negative indices count from the end of the list: 


In [101] : = V [ [-1] ] 
Out [101] = 7+1 


Lists can contain arbitrary elements (for example strings, graphics, expressions, lists, 
functions, etc.). 

If two vectors a and b of equal length are defined, then their scalar product a* ■ b 
is calculated with 


1.7. VECTORS AND MATRICES 


25 


in [ 102 ] : = a = {0.1, 0.2, 0.3 + 2 I>; 
in[io3] :=b = {-0.27 I, 0, 2}; 
in[i04] : = Conjugate [a] ,b 
Out [ 104 ] = 0 . 6 - 4.027 I 


Vectors can be element-wise added, subtracted, multiplied etc. with the usual oper- 
ators: 

In[105] : = a + b 

Out [ 105 ] = {0.1 - 0.27 I, 0.2, 2.3 + 2. 1} 

In[ 106 ] : = 2 a 

Out [106] = {0.2, 0.4, 0.6 + 4. 1} 


1.7.2 matrices 

http://reference.wolfram.com/mathematica/tutorial/BasicMatrixOperations.html 

Matrices are lists of lists, where each sublist describes a row of the matrix: 

in[io7] :=M = {{3, 2, 7}, {1,1, 2}, {0,-1, 5}, {2, 2, 1»; 
in[io8] :=Dimensions [M] 

Out [108] = {4, 31- 


In this example, M is a 4 x 3 matrix. Pretty-printing a matrix is done with the Matrix- 
Form wrapper, 

in [ 109 ] :=MatrixForm[M] 


Accessing matrix elements is analogous to accessing vector elements: 

In[110] : = M [[1,3]] 
out [no] = 7 
In[lll] : = M [ [2] ] 

Out [in] = { 1 , 1, 21- 


Matrices can be transposed with Transpose [M] . 

Matrix-vector and matrix-matrix multiplications are done with the . operator: 

In[112] : = M . a 

0ut[U2] = {2.8 + 14. I, 0.9 + 4. I, 1.3 + 10. I, 0.9 + 2. 1} 


1.7.3 sparse vectors and matrices 

http : / / reference . wolfram . com/language/guide/SparseArrays . html 

Large matrices can take up enormous amounts of computer memory. But in 
practical situations we are often dealing with matrices which are “sparse”, meaning 
that most of their entries are zero. A much more efficient way of storing them is 
therefore as a list of only their nonzero elements, using the SparseArray function. 

A given vector or matrix is converted to sparse representation with 
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In[113] : = M = {{0,3, 0,0, 0,0, 0,0, 0,0}, 
{ 0 , 0 , 0 ,- 1 , 0 , 0 , 0 , 0 , 0 , 0 }, 
{ 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 }}; 
in[ii4] :=Ms = SparseArray [M] 
out [ii4] = Spar seArray [<2> , {3, 10}] 


where the output shows that Ms is a 3 x 10 sparse matrix with 2 non-zero entries. We 
could have entered this matrix more easily by giving the list of non-zero entries, 


in [115] :=Ms = SparseArray [{{1 , 2} -> 3, {2, 4} -> -1}, {3, 10}]; 


which we can find out from 


in [ii6] : = ArrayRules [Ms] 

Out [116] = {{1 , 2} -> 3, {2, 4} -> -1, {_, _} -> 0} 


which includes a specification of the default pattern {_ , _}. This sparse array is con- 
verted back into a normal array with 


in [ii7] :=Normal [Ms] 

Out [117] = {{0,3, 0,0, 0,0, 0,0,0, 0} , 
{ 0 , 0 , 0 ,- 1 , 0 , 0 , 0 , 0 , 0 , 0 }, 
{ 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 }} 


Sparse arrays and vectors can be used just like full arrays and vectors (they are in- 
ternally converted automatically whenever necessary). But for some linear algebra 
operations they can be much more efficient. A matrix multiplication of two sparse 
matrices, for example, scales only with the number of non-zero elements of the ma- 
trices, not with their size. 

1.7.4 matrix diagonalization 

“Solving” the time-independent Schrodinger equation, as we will be doing in sec- 
tion 2.2, involves calculating the eigenvalues and eigenvectors of Hermitian 4 matri- 
ces. 

In what follows it is assumed that we have defined H as a Hermitian matrix. As 
an example we will use 


{{0, 

0.3, 

I, 

0}, 

{0.3 

, 1, 

0 , 

0}, 

{-I, 

0 , 

1. 

-0.2}, 

{0, 

0 , 

-0.2, 

3}}; 


4 A complex matrix H is Hermitian if H = if 4 . See http://en.wikipedia.org/wiki/Hermitian 
matrix. 
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eigenvalues 

The eigenvalues of a matrix H are computed with 
in[ii 9 ] ^Eigenvalues [H] 

out [ii9] = {3.0237, 1.63842, 0.998322, -0.660442} 


Notice that these eigenvalues (energy values) are not necessarily sorted, even though 
in this example they appear in descending order. For a sorted list we use 

in[i 20 ] : = Sort [Eigenvalues [H] ] 

out [i20] = {-0. 660442, 0.998322, 1.63842, 3.0237} 


For very large matrices H, and in particular for sparse matrices (see subsection 1.7.3), 
it is computationally inefficient to calculate all eigenvalues. Further, we are often 
only interested in the lowest-energy eigenvalues and eigenvectors. There are very 
efficient algorithms for calculating extremal eigenvalues, 5 which can be used by 
specifying options to the Eigenvalues function: if we only need the largest two 
eigenvalue, for example, we call 


in[i2i] —Eigenvalues [H, 2, Method 

-> {"Arnoldi", 


"Criteria" -> "RealPart", 


Maxlterations -> 10~6}] 

out [i2i] = {3.0237, 1.63842} 



There is no direct way to calculate the smallest eigenvalues; but since the smallest 
eigenvalues of H are the largest eigenvalues of -H we can use 


in[i 22 ] : = -Eigenvalues [-H, 2, Method 

-> {"Arnoldi", 


"Criteria" -> "RealPart", 


Maxlterations -> 10~6}] 

out [ 122 ] = {0.998322, -0.660442} 



eigenvectors 

The eigenvectors of a matrix H are computed with 


in [ 1 23 ] —Eigenvectors [H] 

out[i23] = {{0. -0.03946131, 0. -0 . 005849891 , -0.117564, 0.992264}, 
{0. +0.5336421, 0 . +0 . 2507621 , 0.799103, 0.117379}, 

{0. -0.00534721, 0. +0 . 9559231 , -0.292115, -0.029187}, 
{0. -0.8447721, 0 . +0 . 1526291 , 0.512134, 0.0279821}} 


In this case of a 4 x 4 matrix, this generates a list of four 4-vectors which are ortho- 
normal. 

Usually we are interested in calculating the eigenvalues and eigenvectors at the 
same time: 

5 Arnoldi-Lanczos algorithm: http://en.wikipedia.org/wiki/Lanczos_algorithm 
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in [ 124 ] :=Eigensystem [H] 

out [124] = {{3.0237, 1.63842, 0.998322, -0.660442}, 

{{0. -0.03946131, 0. -0.005849891, -0.117564, 0.992264}, 
-[0. +0.5336421, 0 . +0 . 2507621 , 0.799103, 0.117379}, 

{0. -0.00534721, 0 . +0. 9559231 , -0.292115, -0.029187}, 
-CO. -0.8447721, 0 . +0 . 1526291 , 0.512134, 0.0279821}}} 


which generates a list containing the eigenvalues and the eigenvectors. The ordering 
of the elements in the eigenvalues list corresponds to the ordering in the eigenvec- 
tors list; but the sorting order is generally undefined. To generate a list of (eigen- 
value, eigenvector) pairs in ascending order of eigenvalues, we calculate 


in [ 125 ] : = Sort [Transpose [Eigensystem[H] ] ] 
Out [ 125 ] = {{— 0 . 660442 , 


{0. -0.8447721, 0 . +0 . 1526291 , 

0.512134, 0.0279821}}, 

-[0.998322, 


{0. -0.00534721, 0. +0.9559231 

-0.292115, -0.029187}}, 

-[1.63842, 


{0. +0.5336421, 0 . +0 . 2507621 , 

0.799103, 0.117379}}, 

{3.0237, 


{0. -0.03946131, 0. -0.005849891, -0.117564, 0.992264}}} 


To generate a sorted list of eigenvalues eval and a corresponding list of eigenvectors 
evec we calculate 


in[i 26 ] : = {eval , evec} = Transpose [Sort [Transpose [Eigensystem[H] ]] ] ; 
In [127] : = eval 

Out [127] = {-0.660442, 0.998322, 1.63842, 3.0237} 

In [128] : = evec 

out [ 128 ] = {{0. -0.8447721, 0 . +0 . 1526291 , 0.512134, 0.0279821}, 

-[0. -0.00534721, 0 . +0 . 9559231 , -0.292115, -0.029187}, 

{0. +0.5336421, 0. +0 . 2507621 , 0.799103, 0.117379}, 

{0. -0.03946131, 0. -0.005849891, -0.117564, 0.992264}} 


The trick with calculating only the lowest- energy eigenvalues can be applied to eigen- 
value calculations as well, since the eigenvectors of -H and H are the same: 


in[i29] : = {eval , evec} = Transpose [Sort [Transpose [-Eigensystem[-H, 2, 
Method -> {"Arnoldi", "Criteria" -> "RealPart", 
Maxlterations -> 10~6}] ] ] ] ; 

in [130] : = eval 

0ut[i30] = {-0. 660442, 0.998322} 

In [131] : = evec 

Out [131] = {{-0 . 733656+0 .4187941 , 0 . 132553-0 . 07566561 , 

-0.253889-0.4447711, -0.0138721-0.0243015 I}, 

-[-0 . 000575666-0 . 005316121 , 0 . 102912+0 . 9503671 , 

-0 . 290417+0 . 03144841 , -0 . 0290174+0 . 00314221}} 


1.8. COMPLEX NUMBERS 


29 


Notice that these eigenvectors are not the same as those calculated further above! 
This difference is due to arbitrary multiplications of the eigenvectors with phase fac- 
tors e" s> . 

To check that the vectors in evec are ortho-normalized, we calculate the matrix 
product 


in[i32] : = Conjugate [evec] . Transpose [evec] // Chop // MatrixForm 


and verify that the matrix of scalar products is indeed equal to the unit matrix. 

To check that the vectors in evec are indeed eigenvectors of H, we calculate all 
matrix elements of H in this basis of eigenvectors: 


in[i33] : = Conjugate [evec] . H. Transpose [evec] // Chop // MatrixForm 


and verify that the result is a diagonal matrix whose diagonal elements are exactly 
the eigenvalues eval. 

1.7.5 exercises 

Q1.20 Calculate the eigenvalues and eigenvectors of the Pauli matrices: 
http : //en . Wikipedia . org/wiki/Pauli_matrices 
Are the eigenvectors ortho-normal? If not, find an ortho-normal set. 

1 .8 complex numbers 

By default all variables in Mathematica are assumed to be complex numbers, unless 
otherwise specified. All mathematical functions can take complex numbers as their 
input, often by analytic continuation. 

The most commonly used functions on complex numbers are Conjugate, Re, 
Im, Abs, and Arg. When applied to numerical arguments they do what we expect: 


in[i34] : = Conjugate [2 + 3*1] 
Out [134] = 2 - 3*1 
In [135] : = lm[0. 7] 

Out [135] = 0 


When applied to variable arguments, however, they fail and frustrate the inexperi- 
enced user: 


in[i36] : = Conjugate [x+I*y] 

out [136] = Conjugate [x] - I *Con jugate [y] 

In[137] : = Im[a] 

0ut[i37] = Im[a] 


This behavior is due to Mathematica not knowing that x, y, and a in these examples 
are real-valued. There are several ways around this, all involving assumptions. The 
first is to use the ComplexExpand function, which assumes that all variables are real : 
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in[i38] : = Conjugate [x+I*y] // ComplexExpand 
Out [138] =x - I*y 

in [139] : = Im [a] // ComplexExpand 
Out [139] = 0 


The second is to use explicit local assumptions, which may be more specific than 
assuming that all variables are real-valued: 

in [140] :=Assuming [Element [x, Reals] && Element [y, Reals], 

Conjugated + I y] // FullSimplify] 
out [140] =x - I*y 

in[i4i] :=Assuming [Element [a, Reals], Im[a]] 

Out [141] = 0 


The third is to use global assumptions (in general, global system variables start with 
the $ sign) : 


in[i42] : = $Assumptions = Element [x, Reals] && Element [y, Reals] && 
Element [a, Reals]; 

in [143] : = Conjugate [x+I*y] // FullSimplify 
Out [143] =x - I*y 

in [144] : = Im [a] // FullSimplify 
Out [144] = 0 


1.9 units 


http : // reference . wolfram . com/mathemat ica/tutor ial/Unit sOverview . html 

Mathematica is capable of dealing with units of measure, as required for physical 
calculations. For example, we can make the assignment 


in [ 145 ] : = s = Quantity["3 m"] ; 


to specify that s should be three meters. A large number of units can be used, as well 
as physical constants: 


in [146] :=kB = Quant ity [" BoltzmannConst ant "] ; 


will define the variable kB to be Boltzmann’s constant. Take note that complicated or 
slightly unusual quantities are evaluated through the online service Wolfram Alpha, 
which means that you need an internet connection in order to evaluate them. 

If you are unsure whether your expression has been interpreted correctly, the full 
internal form 

in [147] :=FullForm[kB] 

out [147] = Quant ity [1 , "BoltzmannConstant"] 


1.9. UNITS 


31 


usually helps. 

In principle, we can use this mechanism to do all the calculations in this lecture 
with units; however, for the sake of generality (as many other computer programs 
cannot deal with units) when we do numerical calculations, we will convert every 
quantity into dimensionless form in what follows. 

In order to eliminate units from a calculation, we must determine a set of units in 
which to express the relevant quantities. This means that every physical quantity x 
is expressed as the product of a unit xq and a dimensionless multiplier x! . The actual 
calculations are performed only with the dimensionless multipliers. For example, a 
length s = 3m can be expressed with the unit .s'o = 1 nr and s' = 3, such that s' sq = s. 
It can equally well be expressed with Sq = 52.9177pm (the Bohr radius) and s' = 
5.669 18 x 10 10 . A smart choice of units can help in implementing a problem. 

As an example we calculate the acceleration of an A380 airplane [m = 560 1) due 
to its jet engines [F - 4 x 311kN). The easiest way is to use Mathematica’s built-in 
unit processing: 

in[i48] :=F = Quant ity [ "4*3 1 1 kN"] ; 

in [ 149 ] :=m = Quantity ["560 t"] ; 

in[i50] : = a = UnitConvert [F/m, "m/s~2"] //N 

out [ 150 ] = 2. 22143 m/s~2 


Now we do the same calculation without Mathematica’s unit processing. Setting 
F = F'Fq, m - m'nio, and a = a! do, Newton’s equation F - ma can be solved for the 
dimensionless acceleration a ': 


C t XN . V. _L . _L J 

m' mo ao 

SI units: With SI units {Fq = IN, mo - 1kg, do = lm/s 2 ), the last term of Equa- 
tion (1.1) is F 0 I [mo ao) = 1, which simplifies the calculation greatly. This is 
the advantage of SI units. 

With F' = 1244000 and m' = 560000 we find d' = F' lm' = 2.22143. Therefore 
we know that the airplane’s acceleration will be d = a! do — 2.221 43 m/s 2 . 

Arbitrary units: If we choose, for example, the units 

• Fq — 1 000 N, the maximum force a human can apply, as the unit of force, 

• m o = 5g, the weight of a hummingbird, as the unit of mass, 

• do — 9.81 m/s 2 , the earth’s gravitational acceleration, as the unit of accel- 
eration, 

the last term Ic = Fq/ { mo do) of Equation (1.1) is computed in Mathematica 
with 

1 in[i5i] :=F0 = Quantity ["1000 N"] ; 

2 in[i52] :=m0 = Quantity ["5 g"l ; 

3 in [ 153 ] : = a0 = Quant ity [" 9 . 81 m/s~2"] ; 

4 in [ 154 ] :=k = F0/ (m0 a0) 

5 out [ 154 ] = 20387. 4 


and we find the airplane’s acceleration with 
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in [155] : = F = Quant ity [ 11 4*31 1 kN"]/FO //N 
Out [155]= 1244. 

in[i56]:=m = Quant ity [" 560 t"]/mO //N 
Out [156] = 1 . 12*10~8 
In[157] : = a = F/m * k 
out [157] = 0 . 226445 


Thus we know that the acceleration is 0.226445g, which is 

In [158] : = a*aO 

0ut[i58] = 2. 22143 m/s~2 


Chapter 2 

quantum mechanics 


In this chapter we connect quantum mechanics to representations that a computer 
can understand. 


2. 1 basis sets and representations 

Quantum- mechanical problems are usually specified in terms of operators and wave- 
functions. The wavefunctions are elements of a Hilbert space; the operators act on 
such vectors. How can these objects be represented on a computer, which only un- 
derstands numbers but not Hilbert spaces? 

In order to find a computer- representable form of these abstract objects, we as- 
sume that we know an ortho-normal 1 basis (| f)} z - of this Hilbert space, with scalar 
product (i\j) — Sij. In section 2.4 we will talk about how to construct such bases. 
For now we make the assumption that this basis is complete, such that^/ UXfi = 1. 
We will see in subsection 2.1.1 how to deal with incomplete basis sets. 

Given any operator si acting on this Hilbert space, we use the completeness re- 
lation twice to find 


si — i ■ si • i = 


Lux*! 


Lb'Xji 


ij 


( 2 . 1 ) 


If we define a numerical matrix A with elements Aij = ( i \si \ j) £ (D we rewrite this as 


J=Z A U l*>0'l- 

ij 


( 2 . 2 ) 


The same can be done with a state vector |i //): using the completeness relation, 


ty) = i- \v) = 


£U'Xz'l 


■\y/) = Y,d \V) l*>. 


(2.3) 


l l 

and defining a numerical vector if/ with elements i/q = (i\y/) e (D the state vector is 

\y/) = Y j y/i\i). (2.4) 

i 

J The following calculations can be extended to situations where the basis is not ortho-normal. For the 
scope of this lecture we are however not interested in this complication. 
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Both the matrix A and the vector if/ are complex-valued objects which can be rep- 
resented in any computer system. Equation (2.2) and Equation (2.4) serve to convert 
between Elilbert-space representations and number-based (matrix/ vector-based) rep- 
resentations. These equations are at the center of what it means to find a computer 
representation of a quantum- mechanical problem. 

2.1.1 incomplete basis sets 

For infinite-dimensional Hilbert spaces we must usually content ourselves with fi- 
nite basis sets which approximate the low-energy physics (or, more generally, the 
physically relevant dynamics) of the problem. In practice this means that an or- 
thonormal basis set may not be complete: 

!>'><*'I = -P (2.5) 

i 

which is the projector onto that subspace of the full Hilbert space that the basis is 
capable of describing. We denote Q = 1 - P as the complement of this projector: Q is 
the projector onto the remainder of the Hilbert space that is left out of this truncated 
description. The equivalent of Equation (2.1) is then 


J=t-J-i = (P + Q)-J-{P+Q) = P-J-P+P-J-q+q-J-P+q-Jq 


= Y^Aij 10 <J I + P-J-Q + Q-J-P + Q-J-Q 

^ J > neglected coupling to (high-energy) part neglected (high-energy) part 

within described subspace 


(2.6) 


In the same way, the equivalent of Equation (2.3) is 


\y/) = t-U ff) = iP+Q)-W = 


U‘> 


+ Q\yA 

neglected (high-energy) part 


within described subspace 


(2.7) 


Since Q is the projector onto the neglected subspace, the component Q \y/) of Equa- 
tion (2.7) is the part of the wavefunction | if/) that is left out of the description in the 
truncated basis. In specific situations we will need to make sure that all terms in- 
volving Q in Equation (2.6) and Equation (2.7) can be safely neglected. 

2.1.2 exercises 

Q2.1 We describe a spin- 1/2 system in the basis containing the two states 
|| = cos | ^ j ll> + e lip sin |^| ||> 

|| {d,(p)) = - e~ 1(p sin(|)lT) + cos(|)|i) (2.8) 

1 . Show that this basis is orthonormal. 

2. Express the states [ f ) and 1 1) as vectors in this basis. 

3. Express the Pauli operators a x , by, d z as matrices in this basis. 
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4. Show that If {Q,ip)) and If {Q,ip)) are eigenvectors oia{Q,ip) - d^sinfS) cos {cp)+ 
a y sin(f)) sin((/j) + <r z cos($). What are the eigenvalues? 

Q2.2 The eigenstate basis for the description of the infinite square well of unit width 
is made up of the ortho -normalized functions 

{x\n) - (p n {x) - \f2sm(nnx) (2.9) 

defined on the interval [0, 1], with n e {1,2,3, .. .}. 

1. Calculate the function PooU>y) = {x\ |tt)(n|] |y>. 

2. In computer-based calculations we limit the basis set to n e {1,2,3,..., n max ( 

for some large value of n max . Using Mathematica, calculate the function 
P„max( x >y) = Ol \n){n\\ | y) (use the Sum function). Make a plot for 

n max = 100 (use the DensityPlot function). 

3. What does the function P represent? 

2.2 time-independent Schrodinger equation 

The time-independent Schrodinger equation is 

jfc\y/) = E\y/). (2.10) 

As in section 2.1 we use a computational basis to express the Hamiltonian operator 
and the wavefunction i// as 


je=Y^ H ij\i) (j\ 
ij 

W=I>/ 10 

With these substitutions the 

L«ij 

. ij 


Y J H iJ Vj\i)=Y, E V'tM 

ij e 


Schrodinger equation becomes 


10 O' I 


E^fc l fc > 


= £ 


i o 


ZHijVk (j\k) \i) = Y J Eyi i \£) 

ijk e 


( 2 . 11 ) 


( 2 . 12 ) 


Multiplying this equation by (m| from the left, and using the orthonormality of the 
basis set, gives 


10 = {m\Y,Ey/( \£) 


E Hi jy/ j (m | 0 = E Eye(m\£) 

= ^mi = ^m£ 

E H m j\f/ j — Ey/ m 


= S m( 


(2.13) 
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In matrix notation this can be written as 


H if/ = Eif /. 


(2.14) 


This is the central equation of this lecture. It is the time-independent Schrodinger 
equation in a form that computers can understand, namely an eigenvalue equation 
in terms of numerical (complex) matrices and vectors. 

Ifyou think that there is no difference between Equation (2.10) and Equation (2.14), 
then I invite you to re-read this section as I consider it extremely important for what 
follows in this course. You can think of Equation (2.10) as an abstract relationship 
between operators and vectors in Hilbert space, while Equation (2.14) is a numerical 
representation of this relationship in a concrete basis set { | f ) } z - . They both contain 
the exact same information (since we converted one to the other in a few lines of 
mathematics) but they are conceptually very different, as one is understandable by 
a computer and the other is not. 


2.2.1 diagonalization 

The matrix form of Equation (2.14) of the Schrodinger equation is an eigenvalue 
equation as you know from linear algebra. Given a matrix of complex numbers H 
we can find the eigenvalues Ei and eigenvectors if/ , using Mathematica’s built-in 
procedures, as described in subsection 1.7.4. 


2.2.2 exercises 

Q2.3 Express the spin- 1/2 Hamiltonian 

Jff = sin(d) cos(<p)<f x + sin($) sin(<p)<ty + cos(-0)<t z (2.15) 

in the basis {[ T> , II)}, and calculate its eigenvalues and eigenvectors. NB: d x ,y,z 
are the Pauli operators. 

2.3 time- dependent Schrodinger equation 

The time-dependent Schrodinger equation is 

ih^-\y/U))=j&(V\y/U)), (2.16) 

at 

where the Hamiltonian J6 can have an explicit time dependence. This differential 
equation has the formal solution 

|i/df)> =‘% , (£ 0 ; t)\yUo)) (2.17) 


in terms of the propagator 

< %'(toH) = l-p f dfi^Hi)- f At\ f dt 2 J&(t!)J&(t 2 ) 
n Jt 0 nr Jt 0 Jt 0 

+ i f Ah f dt 2 f 

n Jt 0 Jt 0 Jt 0 

+ TI [ dt i f dt 2 [ df 3 [ 3 dt 4 J&(t 1 )J?(t 2 )j£ , (t 3 )J&(t 4 ) + ... (2.18) 

n Jt 0 Jt 0 Jt 0 Jt 0 
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which propagates any state from time to to time t. An alternative form is given by 
the Magnus expansion 2 


°UUg t ) = exp 


5^ ol 0 

k = 1 


(2.19) 


with the contributions 


-f 

J t() 


t) — — — I 

to 

n 2 («o = -2 i f 


dt 2 [^f(ti),^f(t 2 )] 


*0 


n3(to; t) = ^3 J t dh d tz j dt3 [J^(t 2 ),^(t3)]] + [J^(t3), [^(t 2 ),^(ti)]]) 


( 2 . 20 ) 


This expansion in terms of different-time commutators is often easier to evaluate 
than Equation (2.18), especially when the contributions vanish for k > fc max (see 
subsection 2.3.3 for the case k max = 1). Even if higher- order contributions do not 
vanish entirely, they (usually) decrease in importance much more rapidly with in- 
creasing Ic than those of Equation (2.18). Also, even if the Magnus expansion is arti- 
ficially truncated (neglecting higher-order terms), the quantum-mechanical evolu- 
tion is still unitary; this is not the case for Equation (2.18). 


2.3.1 time-independent basis 

We again express the wavefunction in terms of the chosen basis, which is assumed 
to be time-independent. This leaves the time dependence in the expansion coeffi- 
cients, 

jfrit) = Y f Ht ] U)\i)(j\ 
ij 

lydfl) = !>/(*) I«>- (2.2i) 

i 

Inserting these expressions into the time-dependent Schrodinger Equation (2.16) 
gives 


i hY^Vdt) It) = 


L Ho v) 1001 

‘j 


!>*:(« I k) =X]H i7 (0^( t) 10. (2.22) 


Multiplying with (£\ from the left: 


or, in matrix notation, 


ihy/( ( f ) = 52 He jit) y/j ( t ) 

i 


ihy/{t) - Hit ) • ij/it). 


(2.23) 


(2.24) 


Since the matrix H{t) is supposedly known, this equation represents a system of 
coupled complex differential equations for the vector iff if), which can be solved on 
a computer. 


^http : / /en. Wikipedia. org/wiki/Magnus_expansion 
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2.3.2 time-dependent basis: interaction picture 

It can be advantageous to use a time-dependent basis. The most frequently used 
such basis is given by the interaction picture of quantum mechanics, where the 
Hamiltonian can be split into a time-independent principal part and a small 
time-dependent part J€\\ 

je{i) = jeQ + je x {f). ( 2 . 25 ) 

Assuming that we can diagonalize possibly numerically, such that the eigen- 
functions satisfy J€§\i) - E;\i), we propose the time -dependent basis 

|i(f)> = e~ iE ‘ tlh \i). (2.26) 

If we express any wavefunction in this basis as 

WU)) = £>/( t) |i(*)> - Y,VMe~ iEitlh \i), (2.27) 

i i 

the time-dependent Schrodinger equation becomes 

£ [i hy/tU) + Etyrdt)] e~ iE ‘ tlh \i) = Y j y/j{t)e~ iE J tlh E j \ j) + Y j y/jU)e~ iE l tlh ^ 1 U) \j) 
i j J 

Y.ihysdUe-^'^i) = i{t)e-' H i tlh .X\(t) \ j) (2.28) 


{k\^VjWe- iE i tlh ^{t)\j) 
j 

Y^VjWe~ iEitlh {k\^iU)\j) 

j 

Y j y/ j U)e~ iiB i~ E ^ tlH {k\J^iU)\j). (2.29) 

j 

This is the same matrix/vector evolution expression as Equation (2.24), except that 
here the Hamiltonian matrix elements must be defined as 

HijU ) = (imt)\j)e- iiE r E ‘» lh . (2.30) 

We see immediately that if the interaction Hamiltonian vanishes (t) = 0], then 
the expansion coefficients i/q(f) become time-independent, as expected since they 
are the coefficients of the eigenfunctions of the time-independent Schrodinger equa- 
tion. 

When a quantum-mechanical system is composed of different parts which have 
vastly different energy scales of their internal evolution then the use of Equa- 
tion (2.30) can have great numerical advantages. It turns out that the relevant in- 
teraction terms Hij(t) in the interaction picture will have relatively slowly evolving 
phases exp [ - i ( Ej - E{) t/H], on a time scale given by relative energy differences and 
not by absolute energies; this makes it possible to numerically solve the coupled dif- 
ferential equations of Equation (2.24) without using an absurdly small time step. 


*' j 

Multiply by (fc| from the left: 

{k\Y j iKy t {t)e- iEitlh \i) = 
i 

Y j ihy/ i U)e~ iEitlh (k\i) = 

=Ski 

ihy/ k U) = 
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2.3.3 special case: [J€{ t ) , t ‘ ') ] = 0 V ( t, t') 

If the Hamiltonian commutes with itself at different times, [Jt( t ) , t') ] = 0 V(f, t'), 

the propagator (2.19) of Equation (2.16) can be simplified to 


<&(t 0 ; t) = exp 


-j: f J?(s)ds 
n Jto 


(2.31) 


and the corresponding solution of Equation (2.24) is 


y/(t) = exp 


-if H{s)ds 
n Jt 0 


■iff (t 0 ). 


(2.32) 


Notice that exponential in this expression has a matrix as its argument: in Mathe- 
matica this matrix exponentiation is done with the Matr ixExp function. 


2.3.4 special case: time-independent Hamiltonian 

In the special (but common) case where the Hamiltonian is time-independent, the 
integral in Equation (2.32) can be evaluated immediately, and the solution is 


y/(t) - exp 


i (t- t 0 ) 
h 


H 


■ip (t 0 ). 


(2.33) 


If we have a specific Hamiltonian matrix H defined, for example the matrix of 
subsection 1.7.4, we can calculate the propagator U(t) — exp(-i Htlh) with 


i in[i59] :=U[t_] = MatrixExp[-I H t] 


where we have used t = ( t - to)/// by expressing time in units of inverse energy (see 
section 1.9). The resulting expression forU [t] will in general be very long. 

2.3.5 exercises 

Q2.4 Demonstrate that the propagator (2.31) gives a wavefunction (2.17) which sat- 
isfies Equation (2.16). 

Q2.5 Calculate the propagator of the Hamiltonian of Q2.3 (page 36). 

2.4 basis construction 

In principle, the choice of basis set {| *">}/ does not influence the way a computer pro- 
gram like Mathematica solves a quantum-mechanical problem. In practice, how- 
ever, we always need a constructive way to find some basis for a given quantum- 
mechanical problem. A basis which takes the system’s Hamiltonian into account 
may give a computationally simpler description; however, in complicated systems 
it is often more important to find any way of constructing a usable basis set than 
finding the perfect one. 
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2.4.1 description of a single degree of freedom 

When we describe a single quantum-mechanical degree of freedom, it is often pos- 
sible to deduce a useful basis set from knowledge of the Hilbert space itself. This 
is what we will be doing in chapter 3 for spin systems, where the well-known Dicke 
basis {|S, M s )} s Ms= _ s turns out to be very useful. 

For more complicated degrees of freedom, we can find inspiration for a basis 
choice from an associated Hamiltonian. Such Hamiltonians describing a single de- 
gree of freedom are often so simple that they can be diagonalized by hand. If this is 
not the case, real-world Hamiltonians can often be decomposed into a “simple” 
part that is time-independent and can be diagonalized easily, and a “difficult" 
part J€\ that usually contains complicated interactions and/or time-dependent terms 
but is of smaller magnitude: 

j£(f) + (2.34) 

A natural choice of basis set is the set of eigenstates of ,77'o , or at least those eigen- 
states below a certain cutoff energy since they will be optimally suited to describe 
the complete low-energy behavior of the degree of freedom in question. This lat- 
ter point is especially important for infinite-dimensional systems (chapter 4), where 
any computer representation will necessarily truncate the dimensionality, as dis- 
cussed in subsection 2.1.1. 


examples of basis sets for single degrees of freedom: 
spin degree of freedom: Dicke states [ S, M$) 

translational degree of freedom: square-well eigenstates, harmonic oscillator eigen- 
states 

rotational degree of freedom: spherical harmonics 
atomic system: hydrogen-like orbitals 

translation-invariant system: periodic plane waves (reciprocal lattice) 


2.4.2 description of coupled degrees of freedom 


A broad range of quantum-mechanical systems of interest are governed by Hamil- 
tonians of the form 


Jeu) = 


N 

L 

\fc=i 




+ "Z^int(l). 


(2.35) 


where N individual degrees of freedom are governed by their individual Hamiltoni- 
ans (f), while their interactions are described by This is a situation we 

will encounter repeatedly as we construct more complicated quantum-mechanical 
problems from simpler parts. A few simple examples are: 


• A set of N interacting particles: the Hamiltonians ,7C fkl describe the individual 
particles, while J6 i nt describes their interactions. 


• A single particle moving in three spatial degrees of freedom: the Hamiltonians 
,77' fx, - v,zi describe the kinetic energy in the three directions, while ,77; n , con- 
tains the potential energy. 
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• A single particle with internal (spin) and external (motional) degrees of free- 
dom which are coupled through a state-dependent potential in 

The existence of individual Hamiltonians assumes that the Hilbert space 
of the complete system has a tensor-product structure 

y = y cu ® V® ® ■ ■ ■ ® V (JV) , (2.36) 

where each Hamiltonian .7C >k} acts only in a single component space, 

Je {k) = (2.37) 

Further, if we are able to construct bases {|z) (fc) }E for all of the component Hilbert 
spaces V {k \ as in subsection 2.4.1, then we can construct a basis for the full Hilbert 
space V by taking all possible tensor products of basis functions: 

I h, h, . • • , in) = I h> (1) ® I h )' [2) ® • • • ® I i N > tJV) - (2.38) 

This basis will have ]~[E , /?/,- elements, which can easily become a very large number 
for composite systems. 


wave vectors (quantum states) 

A product state of the complete system 

\y/) = ® \y/) i2) ® •• • ® \y/) m (2.39) 

can be described in the following way. First, each single-particle wavefunction is 
decomposed in its own basis as in Equation (2.4), 

n k 

= L (2.40) 

h = 1 

Inserting these expansions into Equation (2.39) gives the expansion into the basis 
functions (2.38) of the full system, 



"1 


"2 



l1//> = 

.*i=l 

® 

E ^ife> (2) 
. *2 = 1 

® ® 

E v%\is) m 

Jn= 1 


«i >12 n N 

= £ .£•••£ K In) (2.41) 

i\ = 1 2*2 = 1 1 


In Mathematica, such a wavefunction tensor product can be calculated as fol- 
lows. For example, assume that psil is a vector containing the expansion of |i f/) m 
in its basis, and similarly for psi2 and psi3. The vector psi of expansion coeffi- 
cients of the full wavefunction \y/) = ® ® |i/z} t3) is calculated with 


i in[i60] :=psi = Flatten [KroneckerProduct [psi 1 , psi2, psi3]] 


See Equation (2.45) for a numerical example as an exercise. 
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operators 

If the Hilbert space has the tensor-product structure of Equation (2.36), then the 
operators acting on this full space are often given as tensor products as well, 

A= a {1) ®a {2) ®...® a {N) , (2.42) 

or as a sum over such products. If every single-particle operator is decomposed in 
its own basis as in Equation (2.2), 

n k n k 

= E E i® (2 - 43) 

**=i/'*=i 

inserting these expressions into Equation (2.42) gives the expansion into the basis 
functions (2.38) of the full system, 



n\ n\ 


ri2 112 


A = 

E E *% h \ u m rn m 

L *1=171=1 


e e 

.12 = 1/2 = 1 

<g). . .(g) 


n N n N 


V V a m . | i N ) 

^ t—* In, In' jv/ 


riN n\ n,2 


n N 


E E-.E E E-.E 

*1 = 1 *2=1 *iV — 1 7l = 1 72 = 1 fjv=l 


a™, -a™ 

*1./1 i2j2 lN,]N\ 


, . , tVjN 

tV=l/N =1 

I l'l, *2, • • • , *JV> O’l » J2> • • • > jAtl • 

(2.44) 


In Mathematica, such an operator tensor product can be calculated similarly to 
In [160] above. For example, assume that al is a matrix containing the expansion 
of a [l) in its basis, and similarly for a2 and a3. The matrix A of expansion coefficients 
of the full operator A — a {[ t « a <2> % cP ] is calculated with 


in[i6i] :=A = KroneckerProduct [al , a2, a3] 


Often we need to construct operators which act only on one of the component spaces, 
as in Equation (2.37). For example, the operator which generalizes the component 
Hamiltonians to the full tensor-product Hilbert space is 


in[i 62 ] :=H1 = KroneckerProduct [hi , 

IdentityMatrix [Dimensions [h2] ] , 
IdentityMatrix [Dimensions [h3] ] ] ; 
in[i63] :=H2 = KroneckerProduct [Ident ityMatrix [Dimens ions [hi] ] , 

h2, 

IdentityMatrix [Dimensions [h3] ] ] ; 
in[i64] :=H3 = KroneckerProduct [IdentityMatrix [Dimensions [hi] ] , 

IdentityMatrix [Dimensions [h2] ] , 
h3] ; 

in [165] :=H = HI + H2 + H3; 


where IdentityMatrix [Dimensions [hi] ] generates a unit matrix of size equal to 
that of hi. In this way, the matrices HI, H2, H3 are of equal size and can be added 
together, even if hi, h2, h3 all have different sizes (expressed in Hilbert spaces of 
different dimensions). 
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2.4.3 exercises 

Q2.6 Two particles of mass m are moving in a three-dimensional harmonic poten- 
tial V{r) = \ ma) 2 r 2 with r = \j x 2 + y 2 + z 2 , and interacting via .s'-wave scatter- 
ing Vi,,, = g<5 3 (r i - r 2 ). 

1. Write down the Hamiltonian of this system. 

2. Propose a basis set in which we can describe the quantum mechanics of 
this system. 

3. Calculate the matrix elements of the Hamiltonian in this basis set. 

Q2.7 Calculate psi in In [160] (page 41) without using KroneckerProduct, but 
using the Table command instead. 

Q2.8 Calculate A in In [161] (page 42) without using KroneckerProduct, but using 
the Table command instead. 

Q2.9 Given two spin-1/2 particles in states 

|i/r> (1) = 0.8|f>- 0.6||> 

|i//> t2) = 0.6i|t> + 0.8|[), (2.45) 

use the KroneckerProduct function to calculate the joint state \iy) = (i//} !li ® 
ly/) <2 \ and compare the result to a manual calculation. In which order do the 
coefficients appear in the result of KroneckerProduct? 
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Chapter 3 


spin systems 


In this chapter we put everything we have studied so far together — Mathematica, 
quantum mechanics, computational bases, units — to study quantum-mechanical 
systems with finite-dimensional Hilbert spaces. Spin systems are the simplest kind 
of such systems. 


3.1 quantum-mechanical spin and angular momentum 
operators 

As you know, quantum mechanics is not limited to spins (angular momentum) of 
length S = 1/2. A spin (angular momentum) of length S, with2S e No, is represented 
most easily in the “Dicke basis” of states |S, Ms) with Ms £ {S, S - 1, S - 2, . . . , -S + 
1, -S}. In what follows we will write M instead of Ms whenever no confusion is pos- 
sible. The operators representing such a spin have the properties 


S+\S,M) = VS[S+ 1) -M(M+ 1) |S,M+ 1> 

(3. l)a 

S- 1 S, M) = %/S{S+ 1) -M(M- 1) [S.M- 1> 

(3.1)b 

S Z \S,M) = M \S,M) 

(3. l)c 

S± — Sx i i5 y 

(3.1)d 


In Mathematica we represent these operators in the Dicke basis as follows, with the 
elements of the basis set ordered with decreasing projection quantum number M: 

in[i66] : = SpinQ [ S ] := IntegerQ [2S] && S>=0 

in[i67] : = splus [0] = {{0}-}- // SparseArray; 

In[ 168 ] : = splus [S_?SpinQ] := splus [S] = 

SparseArray [Band [{1 ,2}] -> Table [Sqrt [S(S+1) -M(M+1) ] , 
{M,S-1,-S,-1}] , {2S+1,2S+1}] 
in[i 69 ] : = sminus [S_?SpinQ] := Transpose [splus [S] ] 
in[i 70 ] : = sx [S_?SpinQ] := sx[S] = (splus [S] +sminus [S] ) /2 
in [i7i] : = sy [S_?SpinQ] := sy[S] = (splus [S] -sminus [S] )/ (21) 

In[172] : = sz [S_?SpinQ] := sz[S] = 

SparseArray [Band [{1 , 1}] -> Range [S , -S , -1] , {2S+1 ,2S+1}] 
in[i73] : = SparseIdentityMatrix [n_Integer/ ;n>=l] : = 
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SparseArray [Band [{1 , l}] -> 1, {n,n}] 
in[i74] : = id [S_?SpinQ] := id[S] = SparseldentityMatrix [2S+1] 


• Notice that we have defined all these matrix representations as sparse matri- 
ces (see subsection 1.7.3), which will make larger calculations much more ef- 
ficient later on. 

• The function SpinQ [S] yields True only if S is a nonnegative half-integer 
value and can therefore represent a physically valid spin. In general, functions 
ending in . . ,Q are questions on the character of an argument: IntegerQ, 

PrimeQ, MemberQ, NumericQ, EvenQ, etc. See 

http: //ref erence .wolfram. com/mathematica/tutorial/PuttingConstraintsOnPatterns .html 

for more information. 

• The operator S+ , defined with splus [S] , contains only one off-diagonal band 
of non-zero values. The Spar seArray matrix constructor allows building such 
banded matrices by simply specifying the starting point of the band and a vec- 
tor with the elements of the nonzero band. 

• The operator S z , defined with sz [S] , shows you the ordering of the basis ele- 
ments since it has the projection quantum numbers on the diagonal. 

• The IdentityMatrix function returns a full matrix, which is not suitable for 
large-scale calculations. It is more efficient to define an equivalent SparseldentityMatr 
function which returns a sparse identity matrix of desired size. 

• The last operator id[S] is the unit operator operating on a spin of length S, 
and will be used below for tensor-product definitions. 

• All these matrices can be displayed with, for example, 


i in [ 175 ] :=sx [3/2] // MatrixForm 


3.1.1 exercises 

Q3.1 Verify that for S = 1/2 the above Mathematica definitions give the Pauli matri- 
ces: Sj = \<j i for i = x,y,z. 

Q3.2 Verify in Mathematica that S^. + Sy + S 2 Z = S(S+ 1)1 and [ S x , Sy] = iS z for several 
values of S. What is the largest value of S for which you can do this verification 
within one minute on your computer? Hint: use the Timing function. 

Q3.3 The operators S X: y iZ are the generators of rotations: a rotation by an angle 
a around the axis given by a normalized vector n is done with the operator 
Rji(a) = exp(-i an - S). Set n = {sin(#)cos(<p),sin(fi)sin(<p),cos(-0)} and calcu- 
late the operator Rn(a) explicitly for S - 0, S = 1 12, and S = 1. Check that for 
a = 0 you find the unit operator. 


3.2. SPIN- 1/2 ELECTRON IN A DC MAGNETIC FIELD 
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3.2 spin- 1/2 electron in a dc magnetic field 


As a first example we look at a single spin S — 1/2. As usual we use the basis con- 
taining the two states | f ) = | \ \ ) and | J.) = [ |, - g), which we know to be eigenstates 
of the operators S 2 and S z . The matrix expressions of the operators relevant for this 
system are given by the Pauli matrices divided by two, 








(3.2) 


In Mathematica we enter these as 

in[i76] : = Sx = sx[l/2]; Sy = sy[l/2]; Sz = sz[l/2]; 


using the general definitions of angular momentum operators given in section 3.1. 
Alternatively, we can write 

in[i77] : = -[Sx,Sy,Sz}- = (1/2) * Table [PauliMatrix [i] , -Ci,l,3}]; 


As a Hamiltonian we use the coupling of this electron spin to an external mag- 
netic field, — -fL- B. The magnetic moment of the electron is (l — fi\\g v S in terms 
of its spin S, the Bohr magneton /ig = 9.27400968(20) x 10 -24 J/T, and the electron’s 
g-factor g e = -2. 0023193043622(15). 1 The Hamiltonian is therefore 

& = -/i B ge(S x B x + S y B y + S Z B Z ). (3.3) 

In our chosen matrix representation this Hamiltonian is 

H — — PBge(SxB x + SyBy + S Z B Z ) — — —pBge + ^ v | • (3.4) 

We define the electron’s g-factor with 

in[i78] :=ge = UnitConvert [ "ElectronGFactor " ] 
out [178] = -2 . 00231930436 


3.2.1 time-independent Schrodinger equation 

The time-independent Schrodinger equation for our spin- 1/2 problem is, from Equa- 
tion (2.14), 

■{f/-E\f/ (3.5) 

We remember from section 1.9 that most quantities in Equation (3.5) carry physical 
units, which the computer cannot deal with. Replacing B X: y iZ - BqB' x ^ z and E = 
EqE' gives the dimensionless equation 


:Mb ge 


B z 

Br + i B 


y 


B x -iB 

-B z 


/ RbBq \ 
1 Eo ) 


X 


_ge If B' z 

2 J \B' X + iB'y 


B' x -iB'y\ 

-B’z J 


• If/ = E'lj/ 


(3.6) 


with - |k 1, For concreteness we choose the following units: 

1 Notice that the magnetic moment of the electron is anti-parallel to its spin (g e < 0) . The reason for 
this is the electron’s negative electric charge. When the electron spin is parallel to the magnetic field, the 

electron’s energy is higher than when they are anti-parallel. 
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magnetic field: Bo — 1 G, a common unit for atomic calculations 

energy: Eq - h x 1MHz, where h — 6.62606957 x 10 -34 Js is Planck’s constant. It is 
common to express energies in units of frequency, where the conversion is 
sometimes implicitly done via Planck’s constant. 


We evaluate the numerical prefactor of Equation (3.6) with 


in[i79] :=k = muB*B0/E0 /. 


out [ 179 ] = 1 . 399625 


{muB -> Quantity ["BohrMagnet on"] , 

BO -> Quantity["l Gauss"], 

E0 -> Quantity ["PlanckConstant"] * 
Quant ity["l MHz"]]- 


The fact that this prefactor k comes out to be of order 1 means that we have chosen 
an appropriate set of units. 

We can now define the Hamiltonian in Mathematica, 

in[i 80 ] :=H [Bx_ , By_ , Bz_] = k * (-ge) * (Sx*Bx+Sy*By+Sz*Bz) 


and find its eigenvalues (in units of Eg) and eigenvectors: 


in [i8i] :=Eigensystem [H [Bx,By ,Bz] ] 


As described in subsection 1.7.4 the output is a list with two entries, the first be- 
ing a list of eigenvalues and the second a list of associated eigenvectors. As long as 
the Hamiltonian matrix was hermitian, the eigenvalues will all be real-valued; but 
the eigenvectors can be complex. Since the Hilbert space of this spin problem has 
dimension 2, and the basis contains two vectors, there are necessarily two eigen- 
values and two associated eigenvectors of length 2. The eigenvalues can be called 
E± = +ipBf>ell-B||, or » in our dimensionless formulation, E' ± = +ky ||B'||. The list of 
eigenvalues is given in the Mathematica output as {E'_,E' + }. Notice that these eigen- 
values only depend on the magnitude of the magnetic held, and not on its direction. 
This is to be expected: the choice of the basis as the eigenstates of the S z operator 
was entirely arbitrary, and therefore the energy eigenvalues cannot depend on the 
orientation of the magnetic Held with respect to this quantization axis. Since there 
is no preferred axis in this system, there cannot be any directional dependence. 

The associated eigenvectors are 


V*± = { 


B Z ±||B|| 
B x + i By 


1 }, 


(3.7) 


which Mathematica returns as a list of lists, Notice that these eigenvectors 

are not normalized. 


3.2.2 exercises 

Q3.4 Calculate the eigenvalues (in units of J) and eigenvectors (ortho-normalized) 
of an electron spin in a magnetic Held of 1 T in the x-direction. 

Q3.5 SetB = B[e x sin(H)cos(</?) + eySin(H)sin(<p) + e-cos(H)] and calculate the eigen- 
values and normalized eigenvectors of the electron spin Hamiltonian. 
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3.3 coupled spin systems: 87 Rb hyperfine structure 

Ground-state Rubidium-87 atoms consist of a nucleus with spin I = 3/2, a single va- 
lence electron (spin S — 1/2, orbital angular momentum L- 0, and therefore total 
spin / = 1/2), and 36 core electrons which do not contribute any angular momen- 
tum. In a magnetic field along the z-axis, the effective Hamiltonian of this system 

is- 

& = + hA hfs I- J- p B 5 z (g/l z + g s S z + giL z ), (3.8) 

where h is Planck’s constant, /ig is the Bohr magneton, /Ihfs = 3.417341305452 145(45) GHz 
is the spin-spin coupling constant in the ground state of 87 Rb, gj - +0.000995 1414(10) 
is the nuclear g-factor, g$ = -2.0023193043622(15) is the electron spin g-factor, 
and gi = -0.99999369 is the electron orbital g-factor. 

The first part of Equation (3.8) contains all electrostatic interactions, core 
electrons, nuclear interactions etc. We will assume that the system is in the ground 
state of which means that the electron is in the 5 2 Si/ 2 state. This ground state is 
eight-fold degenerate and consists of the four magnetic sublevels of the I - 3/2 nu- 
clear spin, the two sublevels of the S — 1/2 electronic spin, and the single sublevel of 
the L- 0 angular momentum. The basis for the description of this atom is therefore 
the tensor product basis of a spin-3/2, a spin-1/2, and a spin-0. 

The spin operators acting on this composite system are defined as in subsec- 
tion 2.4.2. For example, the nuclear-spin operator I x is extended to the compos- 
ite system by acting trivially on the electron spin and orbital angular momenta, 

I x i — ► /*®1®1. The electron-spin operators are defined accordingly, for example 
S x <-* 1 ® S x ® 1 ■ The electron orbital angular momentum operators are, for example, 

L x !-► 1 ® 1 ® L x . In Mathematica these operators are defined with 


In [182] 

= Ix 

= KroneckerProduct [sx [3/2] , 

id[l/2] , 

id[0] ] 


In [183] 

-iy 

= KroneckerProduct [sy [3/2] , 

id[l/2] , 

id[0] ] 


In [184] 

= Iz 

= KroneckerProduct [sz [3/2] , 

id[l/2] , 

id[0] ] 


In [185] 

= Sx 

= KroneckerProduct [id [3/2] , 

sx [1/2] , 

id[0] ] 


In [186] 

= Sy 

= KroneckerProduct [id [3/2] , 

sy [1/2] , 

id[0] ] 


In [187] 

= Sz 

= KroneckerProduct [id [3/2] , 

sz [1/2] , 

id[0] ] 


In [188] 

= Lx 

= KroneckerProduct [id [3/2] , 

id[l/2] , 

1 — 1 
1 1 

0 

1 1 

X 

w 


In [189] 

=Ly 

= KroneckerProduct [id [3/2] , 

id[l/2] , 

sy [0] ] 


In [190] 

= Lz 

= KroneckerProduct [id [3/2] , 

id[l/2] , 

l — l 
l — l 

0 

1 1 

N 

W 



The total electron angular momentum is J - S+L: 


in[i9i] : = Jx = Sx + Lx ; Jy = Sy + Ly; Jz = Sz + Lz; 


The total angular momentum of the 8 ' Rb atom is F - I + J: 


in[i92] :=Fx = Ix + Jx; Fy = Iy + Jy; Fz = Iz + Jz; 


From these we can define the hyperfine Hamiltonian with magnetic field in the 
z-direction as 


2 see http : //steck .us/alkalidata/rubidium87numbers . pdf 
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in [193] :=Hhf = A (lx . Jx+Iy . Jy+Iz . Jz) - muB*Bz*(gI*Iz+gS*Sz+gL*Lz) ; 
in [194] : =hf c = {A -> 3417.341305452145, 
gS -> -2.0023193043622, 
gL -> -0.99999369, 
gl -> +0.0009951414, 
muB -> 1.3996255481168427}; 


where we have made the following assumptions: 

• Energies are expressed in units of MHz, after dividing by Planck’s constant; 
magnetic field strengths are expressed in units of Gauss. This is an alterna- 
tive description of what we did with the constant k in subsection 3.2.1: essen- 
tially we choose a compatible system of units which gives k = 1 (just like the SI 
units) . 

• A = 7l h f s /MHz = 3417.34 

• muB = pv,lh x G/MHz = 1.399625: 

1 in [ 195 ] :=UnitConvert ["BohrMagneton/PlanckConstant" , "MHz/G"] 

2 out [ 195 ]= 1 . 399625 MHz/G 


This yields the Hamiltonian as an 8 x 8 matrix, and we can calculate its eigenvalues 
and eigenvectors with 

in [196] :={eval , evec} = Eigensystem[Hhf ] // FullSimplify ; 


We plot the energy eigenvalues with 

in[i97] :=Plot [Evaluate [eval /. hfc] , {Bz, 0, 3000}, 

Frame -> True, FrameLabel -> {"Bz / G", "E / MHz"}] 



B Z IG 


3.3.1 eigenstate analysis 

In this section we analyze the results eval and evec from the Hamiltonian diagonal - 
ization above. For this we first need to define ortho-normalized eigenvectors since 
in general we cannot assume evec to be ortho-normalized. 

In general we can always define an ortho-normalized eigenvector set with 
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in[i98] :=nevec = Orthogonalize [evec] 


The problem with this definition is, however, immediately apparent if you look at 
the output given by Mathematica: since no assumptions on the reality of the vari- 
ables were made, the orthogonalization is done in too much generality and quickly 
becomes unwieldy. Even using Assuming and ComplexExpand, as in section 1.8, 
does not give satisfactory results. But if we notice that the eigenvectors in evec are 
all purely real-values, and are already orthogonal, then a simple vector-by-vector 
normalization is sufficient for calculating an ortho -normalized eigenvector set: 


in[i99] :=nevec = #/Sqrt[#.#] & /<§ evec; 

in[200] :=nevec . Transpose [nevec] // FullSimplify 


The fact that In [200] finds a unit matrix implies that the vectors in nevec are ortho- 
normal. 


field-free limit 

In the field-free limit B z - 0 the energy levels are 


in[20i] :=Assuming [A > 0, Limit [eval, Bz -> 0]] 

out[20i] = -[3A/4, 3A/4 , -5A/4, 3A/4, -5A/4, 3A/4, -5A/4, 3A/4> 


We see that the level with energy - 1 A is three-fold degenerate while the level with 
energy | A is five-fold degenerate. This is also visible in the eigenvalue plot above. 
Considering that we have coupled two spins of lengths I - 8 and J - |, we expect 
the composite system to have either total spin F - 1 (three sublevels) or F = 2 (five 
sublevels) ; we can make the tentative assignment that the F - 1 level is at energy 
Ei — - 1 A and the F — 2 level at E-i - | A. 

In order to demonstrate this assignment we express the matrix elements of the 
operators F 2 and F z in the field-free eigenstates, making sure to normalize these 
eigenstates before taking the limit B z —> • 0: 


in[202] —nevecO = Assuming[A > 0, Limit [nevec , Bz -> 0] ] ; 
in[203] :=nevec0 . (Fx.Fx+Fy.Fy+Fz.Fz) . Transpose [nevecO] 
in[204] :=nevec0 . Fz . Transpose [nevecO] 


Notice that in this calculations we have used the fact that all eigenvectors are real, 
which may not always be the case for other Hamiltonians. We see that the field- 
free normalized eigenvectors nevecO are eigenvectors of both F 2 and F z , and from 
looking at the eigenvalues we can identify them as 


{|2,2>, |2, -2>, |1,0), |2,0),|1, 1>, [2, 1), |1,-1>, |2,-1>] (3.9) 

in the notation | F, Mp) . These labels are often used to identify the energy eigenstates 
even for B z ^ 0. 
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low-field limit 

For small magnetic fields, we series-expand the energy eigenvalues to first order in 

B z : 

in[205] :=Assuming [A > 0, Series [eval, {Bz, 0, 1}] // FullSimplify] 


From these low-field terms, in combination with the field-free level assignment, we 
see that the F = 1 and F - 2 levels have effective g-factors of gi = - (gs - 5 gj) 14 ~ 
-0.501824 and g 2 = (gs + 3g;)/4 ~ 0.499833, respectively, so that their energy eigen- 
values follow the form 


E f ,m f (B z ) - Ep{0)- (3.10) 


These energy shifts due to the magnetic field are called Zeeman shifts. 

high-field limit 

The energy eigenvalues in the high-field limit are infinite; but we can calculate their 
lowest-order series expansions with 

in[ 206 ] :=Assuming [muB > 0 && gS < -gi < 0, 

Series [eval, {Bz, Infinity, 0}] // FullSimplify] 


From these expansions we can already identify the states in the eigenvalue plot 
above. 

In order to calculate the eigenstates in the high-field limit we must again make 
sure to normalize the states before taking the limit B z — > • oo: 

in[207] := neve c inf = Assuming [A > 0 && muB > 0 && gS < -gi < 0, 

Limit [nevec, Bz -> Infinity]] 


Out [207] = {{ 1 , 

0, 

0, 

0, 

0, 

0, 

0, 

0>, 

{0, 

0, 

0, 

0, 

0, 

0, 

0, 

1>. 

{0, 

0, 

0, 

-1. 

0, 

0, 

0, 

0>, 

{0, 

0, 

0, 

0, 

1. 

0, 

0, 

0>, 

{0, 

-1 

, 0, 

0 

, 0, 

0, 

0, 

0>, 

{0, 

0, 

1. 

0, 

0, 

0, 

0, 

0>, 

{0, 

0, 1 

0, 

0, ' 

0, 

-1. 

0, 

0>, 

{0, 

0, 

0, 

0, 

0, 

0, 

1, 

o» 


From this we immediately identify the high-field eigenstates as our basis functions 
in a different order, 




2 1 2 


2 } 2 ' 


2 } 2 ' 


(3.11) 


where we have used the abbreviation | Mi, Mj) = 1 1, Mj> <s> | Mj). You can verify this 
assignment by looking at the matrix elements of the I z and J z operators with 


in [ 208 ] : = neve c inf 

Iz 

Transpose [nevecinf ] 

in [ 209 ] : = neve c inf 

Jz 

Transpose [nevecinf] 
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3.3.2 “magic” magnetic field 

The energy eigenvalues of the low-field states 1 1, - 1) and |2, 1} have almost the same 
first-order magnetic field dependence since gi ~ -g 2 (see low-field limit above). If 
we plot their energy difference as a function of magnetic field we find an extremal 
point: 

in[ 2 io] :=Plot [eval [ [6] ] -eval [ [7] ] -2A /. hfc, {Bz, 0, 6}] 



B z l G 

At the “magic” field strength Bq - 3.228 95 G the energy difference is independent of 
the magnetic field (to first order): 


in [2 ii] : =F indMinimum [eval [ [6] ] -eval [ [7] ] -2A /. hfc, {Bz, 2}] 
out [ 211 ] = {-0.00449737, {Bz -> 3.22895}} 


3.3.3 coupling to an oscillating magnetic field 

In this section we briefly study the coupling of a 87 Rb atom to a weak oscillating 
magnetic field. Such a field could be the magnetic part of an electromagnetic wave, 
whose electric field does not couple to our atom in the electronic ground state. This 
calculation is a template for more general situations where a quantum-mechanical 
system is driven by an oscillating field. 

The 87 Rb hyperfine Hamiltonian in the presence of an oscillating magnetic field 
is 

jfc(t) = JiA hk I- J- pnB-igJz + gsSz + g L L z ) - cos (cot) X p B B dt • ( gli + gsS + g L L ) 
v , ' „ ' 

(3.12) 

where the static magnetic field is assumed to be in the z direction, as before. Unfor- 
tunately, [^(t),^(t')] = {je x ,jed (cos (cot) - cosfiuf')) ^ 0 in general, so we cannot 
use the exact solution of Equation (2.32) of the time-dependent Schrodinger equa- 
tion. In fact, the time-dependent Schrodinger equation of this system has no ana- 
lytic solution at all. In what follows we will calculate approximate solutions. 

Since we have diagonalized the time-independent Hamiltonian . : 7C {) already, we 
use its eigenstates as a basis for calculating the effect of the oscillating perturbation 
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(f). In general, calling { | / > } / the set of eigenstates of 3€o> with Mo\i) = Ei\i), we 
expand the general hyperfine state as 

mt)) = Y J Vdt)e- iBitlh \i). (3.13) 

i 

The time-dependent Schrodinger equation for the expansion coefficients i/q(t) in 
this interaction picture is given in Equation (2.29): 


ihyrdt) = ]T y/j cos{cot)(i\ J€\ \ j ) 

i 


z j 


-il^E <*>\t i| 

+ e 




Tij, (3.14) 


where we have replaced cos(wf) = \e lwt + \e '" >l and defined 


Tij = {imj) = -{i\ p B B Agil + gsS + gMUj). (3.15) 


From Equation (3.14) we make two key observations: 

Transition matrix elements: The time-independent matrix elements Tij of the per- 
turbation Hamiltonian are called the transition matrix elements and describe 
how the populations of the different eigenstates of J€o are coupled through 
the oscillating field. We calculate them in Mathematica as follows: 


2 

3 

4 

5 

6 

7 

8 
9 

10 


in [ 212 ] : = HO = A (Ix.Jx + Iy.Jy + Iz.Jz) 

- muB Bz (gS Sz + gL Lz + gl Iz) ; 
in[2i3] :=H1 = -muB (gS (Bacx Sx + Bacy Sy + Baez Sz) 

+ gl (Bacx lx + Bacy Iy + Baez Iz) 

+ gL (Bacx Lx + Bacy Ly + Baez Lz) ) ; 

in[2i4] :=H [t_] = HO + HI Cos [w t] ; 

in[2i5] : = {eval , evec} = Eigensystem[HO] // FullSimplify ; 
in[ 2 i 6 ] : = nevec = Map [#/Sqrt [# . #] &, evec]; 
in[2i7] :=T = Assuming [A > 0, 

nevec. HI .Transpose [nevec] // FullSimplify]; 


Looking at this matrix T we see that not all energy levels are directly coupled 
by an oscillating magnetic field. For example, T \ : 2 = 0 indicates that the pop- 
ulations of the states 1 1) and [2) can only be indirectly coupled through other 
states, but not directly. 

Numerical solution: We will use the time unit to - 1 ps. Since our unit of energy is 
Eq = h x 1 MHz, the reduced Planck constant takes on the value h-hl (Eq to) = 
hl[hx 1 MHz x 1 ps) = Hlh= 1 / (2 n) . It is important not to forget this factor in 
the time-dependent Schrodinger equation. 

Equation (3.14) is a series of linear coupled differential equations, which we 
can write down explicitly in Mathematica with 

1 in [ 218 ] := deqs = Table [I*/i*Subscript [psi , i] ’ [t] == 

2 1/2 Sum[Subscript [psi , j] [t] *T [ [i , j] ] * 

3 (E~ (-1* ( (eval [ [j] ] -eval [ [i] ] )/?i-w) t) + 
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E~ (I* ( (eval [ [i] ] -eval [ [ j ] ] )/h- w)t) ) , 
{j,8}], -Ci , 8>] /. h -> 1/ (2*Pi) ; 


where w = co to is the frequency of the magnetic field in units of ps -1 . Assuming 
concrete conditions, for example the initial state \y/[t = 0)> = \F — 2, Mp - -2) 
which is the second eigenstate neve c [ [2] ] [see Equation (3.9)], and magnetic 
fields B z = 3.228 95 G, = lmG, By C = Bf ' = 0, and an ac field frequency of 

oj = 2tt x 6.828GHz, we can find the time- dependent state i//(f)} with 


1 in [ 219 ] : = S = NDSolve [Jo in [deqs/ . hf c/ . {Bz->3 . 22895 , Bacx->0 . 001 , 

2 Bacy->0,Bacz->0,w->2*Pi*6828}, 

3 {Subscript [psi , 1] [0] ==0 , Subscript [psi , 2] [0] ==1 , 

4 Subscript [psi , 3] [0] ==0 , Subscript [psi , 4] [0] ==0 , 

5 Subscript [psi , 5] [0] ==0 , Subscript [psi , 6 ] [0] ==0 , 

e Subscript [psi ,7] [0] ==0 , Subscript [psi , 8] [0] ==0}] , 

7 Table [Subscript [psi , i] [t] ,{i , 8 }] , {t, 0, 30}, 

s MaxStepSize->10~ (-5) , MaxSteps->10~7] 


Notice that the maximum step size in this numerical solution is very small 
(10 _5 fo = lOps), since it needs to capture the fast oscillations of more than 
6.8 GHz. As a result, a large number of numerical steps is required, which 
makes this way of studying the evolution very difficult in practice. 

We can plot the resulting populations with 


1 in[ 220 ] :=Plot [Evaluate [Abs [Subscript [psi , 2] [t] /. S[[l]]]~2], 

2 {t, 0, 30}] 



1 in[ 22 i] :=Plot [Evaluate [Abs [Subscript [psi , 7] [t] /. S[[l]]]“2], 

2 {t, 0, 30}] 
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II pS 


We see that the population is slowly and slightly sloshing between -eigenstates 

1 2) « \F= 2,Mf = -2) and |7) « |F = 1,M P = -1> [see Equation (3.9)]. 

Rotating-wave approximation: The time-dependent prefactor exp | -i f E ’ - w| f J + 

[/ . \ I . jj. 

i ' h 1 - oil t\ ofEquation (3.14) oscillates very rapidly unless either 1 h ' - 

Ei~Ej 

WKOor ' h 1 -co ~ 0, where one of its terms changes slowly in time. The 
rotating-wave approximation (RWA) consists of neglecting all rapidly rotating 
terms in Equation (3.14). Assume that there is a single pair of states |i> and 
| j) such that Ej - Ej ~ hco, with Ej > Ej, while all other states have an energy 
difference far from tuo. The RWA thus consists of simplifying Equation (3.14) 
to 

1 i[^J—co}t 

i hy/i(t)~-y/j(t)e\ I Tjj 

1 -if ^7 P—At 

i hif/j(t)~-y/ i [t)e \ I Tp 

ihcj/kit) a 0 for k€ {i,j} (3.16) 


with Tji - T*. This approximate system of differential equations has the exact 

J IJ 

solution 


i VM = e 2 Af 




Tit 


Vm cos ( — j - i 


V jit) = 

Vklt) = V'fc(O) for kt{i,j] 


2 

Tit 


A 


T, 


y/ilO) cos I ^ I + i I - y/i (0) - -r^ y/j (0) | sin I ^ 


Tl T * ' hTl 
A T h 

-yrjl 0 ) + ^-( 0 ) 


Tit 


2 

Tit 

~2 


(3.17) 


in terms of the detuning A = oi — ( Ej - E j ) / h and the generalized Rabi frequency 
Tl = sJ\Tij\ 2 lh 2 + A 2 . We can see that the population sloshes back and forth 
(“Rabi oscillation”) between the two levels |i) and | j) with angular frequency 
O. 

We can verify this solution im Mathematica as follows. First we define 

3 The following derivation is readily extended to situations where several pairs of states have an en- 
ergy difference approximately equal to hoj. In such a case we need to solve a larger system of coupled 
differential equations. 
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1 in[222] :=Delta = w - (Ei-Ej)/fi; 

2 in [223] := Omega = Sqrt [Tij *Tj i/fi~2 + Delta~2] ; 


and the solutions 

1 in[224] :=psii [t_] = E~ ( - 1 *Delt a*t /2) * (psiiO*Cos [0mega*t/2] 

2 +1* (Delta/Omega*psiiO-Ti j / (fi*0mega) *psi j 0) 

3 *Sin [0mega*t/2] ) ; 

4 in[225] :=psij [t_] = E~ (I*Delta*t/2) * (psij 0*Cos [0mega*t/2] 

5 -I* (Delta/Omega*psi j O+Tj i/ (fi*0mega) *psiiO) 

6 *Sin [0mega*t/2] ) ; 


With these definitions, we can check the Schrodinger equations (3.16): 

1 in[226] :=FullSimplify [I*fi*psii ’ [t] == 

2 (1/2) * psij [t] * Exp[I*((Ei-Ej)/fi-w)*t]*Tij] 

3 out [226]= True 

4 in[227] :=FullSimplify [I*fi*psij ’ [t] == 

s (1/2) * psii [t] * Exp [-1* ( (Ei-Ej ) /fi-w) *t] *Ti j] 

6 out [227]= True 


as well as the initial conditions 

1 In[228] : = psii [0] 

2 0ut[228]=psii0 

3 In[229] : = psij [0] 

4 Out [229] =psi j 0 


dressed states: If we insert the RWA solutions, Equation (3.17), into the definition 
of the general hyperfine state, Equation (3.13), and set all coefficients y/^ = 0 
for k £ {/, j], and then write sin(z) = ( e lz - e _lz )/(2i) and cos(z) = (e lz + e _lz )/2, 
we find the state 


Wit)) ~ y/dt)e lEit,n \ i) + y/jU)e 


- -e 
2 


-i E<- 


B(n-A) 


tin 


l W 


V/(0)|1 + -|-V;(0)^ 


1 


tih 


ho. 


Vh(0)|i--| + ^-(0)^ 


fin 


U) + 

i) + 


A 


j 1 * 

ij 


t,,(°)|l--j-^(°)^ 

^(0)|l+|] + V/(0)^g 


fin 


e lwt \ j) 


e iMt \j) 


(3 


This state consists of two components, called dressed states, which in a more 
complete description (including the quantization of the electromagnetic field 
at angular frequency to ) couple the upper state with n photons in the driving 
field | i, n) and the lower state with n+ 1 photons in the driving field \j,n+ 1). 
The energy of the photon field contributes an additional time dependence 

\i)~\i,n)e~ inu>t , \j) ~ \j,n+ l)e~ i{n+1)u>t , (3.19) 
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and Equation (3.18) becomes 


1 — i^Bi+nH(i)+ 0) j t/h J 

2 6 " | 





^( 0) g]u-, 

-^Sl 1 


n) + 
i, n) + 


A 


— -J — ^(0)^ 


T*. 

‘] 


hn 


I j,n + I) 


A 


V;(0)|l+-| + ^(0)^ 


u 


ho, 


\j.n+ l>Jj. 

(3.20) 


With this substitution, we can see that the dressed states are time-invariant 
apart from their energy prefactors, which correspond to their effective energy 
in the presence of the oscillating field, 


E+ - Ei + nhcj + 


/UA + Q) 
2 


= Ej + {n + 1 )hd) + 


H (- A + Q) 
2 


(3.21) 


We look at these dressed states in two limits: 

• On resonance (A = 0), Equation (3.20) becomes 


|l//(f)> SB 

I e -i(Ei+nHb)~ 5 1 Ttj |) tlh 


+ 1 e -i[Ei+nhu>+ \ | Tij |) t lh 
2 


1/q (0)- 1/^(0) ■ 


Ttj I 


t//i(0) + ^;(0)- 


|i, n) + 
i, n) + 




r j- 

i] 


V;(0) + l//;(0) 


‘‘ J 1 

ij 


\j,n+\) 


Ij.h + D 

(3.22) 


The dressed states are therefore at energies 

£+ = Ei + nhu) ± - 1 Ti i\ = E j + (n + l)hw ± - \T{ ; | (3.23) 

2 2 

in the presence of a resonant ac coupling field. We can say that each en- 
ergy level is split in two, separated by the magnitude \Tij\ of the coupling 
matrix element. The states are equal mixtures of the original states \ i, n) 
and \j,n+l). 

| J 1 . . |2 

• Far off-resonance (A — ► ±oo) we have Q ~ |A|+ 2fe ' 2 J ^ , and Equation (3.20) 
becomes 


[yd?)> ! 


-iiEi+nha)- 


IQjT 

AhA 


tlh 


+ e 


V t (0) \i,ri) 

if hj+(n+\)hco-t \tlh 


yjj(0)\j,n+ 1>. (3.24) 


(Hint: to verify this, look at the cases A — *• +00 and A — * -00 separately). 

| y . .| 2 

The energy levels [ i, n) and \j, n + 1) are thus shifted by + 4f { A , respec- 
tively, and there is no population transfer between the levels. That is, the 
dressed states become equal to the original states. Remember that we 
had assumed Ei > Ef. 
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- For a blue-detuned drive (A — > • + 00 ), the upper level | i) is lowered in 

I Ti -| 2 

energy by A E = while the lower level | j) is raised in energy by 
A E. 

- For a red-detuned drive (A — > • - 00 ), the upper level |/> is raised in 

I T -| 2 

energy by A E - while the lower level \j) is lowered in energy by 

A E. 

These shifts are called ac Zeeman shifts in this case, or level shifts more 
generally. When the oscillating held is a light held, level shifts are often 
called light shifts or ac Stark shifts. 

3.3.4 exercises 

Q3.6 Take two angular momenta, for example 1-3 and J - 5, and calculate the 
eigenvalues of the operators I 2 , I z , J 2 , J z , F 2 , and F z , where F -I + J. 

Q3.7 In Q3.6 you have coupled two angular momenta but you have not used any 
Clebsch-Gordan coefficients. Why not? Where do these coefficients appear? 

Q3.8 For a spin of a certain length, for example S - 100, take the state |S, S) and 
calculate the expectation values (S x ), (S y ), ( S z >, (S 2 )-(S X ) 2 , ( S 2 )-(S y } 2 , (S 2 )- 

<s z ) 2 . 

Q3.9 Show that the results of the numerical solution plotted with In [220] and 
In [221] (page 55) can be reproduced with the RWA solution of Equation (3.17) 
with i - 2 and j — 7. 

Q3.10 Plot the 87 Rb level shifts at B z = 3.228 95 G (the magic held) for the following 
directions of the oscillating magnetic held: 


• circularly polarized around the quantization axis: B‘ k - B{e x + ie v ) 

• linearly polarized parallel to the quantization axis: B'" = Be z 

Which polarizations can be absorbed by 87 Rb at which frequencies? 

Q3.ll Do the presented alkali atom calculation for 23 Na: are there any magic held 
values? 

http : //steck . us/alkalidata/sodiumnumbers . pdf 

Q3.12 Do the presented alkali atom calculation for 85 Rb: are there any magic held 
values? 

http : //steck . us/alkalidata/rubidium85numbers . pdf 

Q3.13 Do the presented alkali atom calculation for 133 Cs: are there any magic held 
values? 

http : //steck . us/alkalidata/cesiumnumbers . pdf 
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3.4 coupled spin systems: transverse Ising model 

We now turn to larger numbers of coupled quantum-mechanical spins. A large class 
of such coupled spin systems can be described with the Hamiltonian 

N N - 1 N 

,#=E^®+E E (3.25) 

k-1 k= 1 k'-k+l 

where the .7C {k> are single-spin Hamiltonians (for example couplings to a magnetic 
field) and the ■% > ^ k ' are coupling Hamiltonians between two spins. Direct cou- 
plings between three or more spins can usually be neglected. 

In particular we study the “transverse Ising" Hamiltonian 

u N N 

Je= — E S® - E S® S? +1) (3-26) 

2 k = 1 k = 1 

acting on a ring of N spin-S systems where the (N + l)st spin is identified with the 
first spin. We can read off three limits from this Hamiltonian: 

• For b — » ±oo the spin-spin coupling Hamiltonian can be neglected, and the 
ground state will have all spins aligned with the ±x direction, 

l^ + oo> = ITx)®E [V-oo> = IL->® W (3.27) 


The system is therefore in a product state for b — *• oo, which means that there 
is no entanglement between spins. In the basis of \S,M) Dicke states, Equa- 
tion (3.1), the single-spin states making up these product states are 


x) - 2 S E 


^ \ 

M=—S \ 


2 S 

M+S 


I S,M), IU=2- S E (-D 
M=-S 


M+S 


N 


2S 

M+S 




\S,M), 

(3.28) 


which are aligned with the x-axis in the sense that S r |] x ) = S|l x ) and ,S', || X ) = 

-SIU>. 


• For b — 0 the Hamiltonian contains only nearest-neighbor ferromagnetic spin- 
spin couplings - S^ 1 S ( t 1 1 ■ We know that this Hamiltonian has two degener- 
ate ground states: all spins pointing up or all spins pointing down, 

l^ 0 T> = ID®E IVoi> = ll z)* N , 0.29) 


where in the Dicke-state representation of Equation (3.1) we have || z ) = |S,+S) 
and I iz) = IS, -S). While these two states are product states, for \b\ « 1 the per- 
turbing Hamiltonian -f XEi S? is diagonal in the states — which 

are not product states. The exact ground state for 0 < b « 1 is close to 

and for -1 « b < 0 it is close to These are both maximally entan- 

gled states (“Schrodinger cat states”). 

Here we calculate the ground-state wavefunction | y/b) as a function of the parameter 
b, and compare the results to the above asymptotic limits. 
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3.4.1 basis set 

The natural basis set for describing a set of N coupled spins is the tensor-product 
basis (see subsection 2.4.2). In this basis, the spin operators S z fc j, jZ acting only on 
spin k are defined as having a trivial action on all other spins, for example 

I— 1 ® 1 ® ■ ■ ■ ® 1 ®S X 8 l8-’'8 1 . (3.30) 

(fc-1) (JV-fc) 


In Mathematica such single-spin-S operators acting on spin k out of a set of N spins 
are defined with 


in[230] : = op [S_?SpinQ , n_Integer, k_Integer, 

a_?MatrixQ] /; 

l<=k<=n && Dimensions [a] == {2S+1 , 2S+1} := 

KroneckerProduct [Sparseldent ityMatrix [ (2S+1 ) ~ (k- 1) ] , 

a J 

Sparseldent ityMatrix [ (2S+1 ) ~ (n-k) ] ] 

in[23i] : = sx [S_?SpinQ , n_Integer, k_Integer] 

/; l<=k<=n := 

op[S, n, k, sx[S]] 


in[232] : = sy [S_?SpinQ , n_Integer, k_Integer] 

/; l<=k<=n := 

op [S , n, k, sy [S] ] 


in[233] : = sz [S_?SpinQ , n_Integer, k_Integer] 

/; l<=k<=n := 

op[S, n, k, sz[S]] 



Notice that we have used n-N because the symbol N is already used internally in 
Mathematica. 

From these we assemble the Hamiltonian: 

in[234] :=H [S_?SpinQ , n_Integer/ ;n>=3 , b_] : = 

-b/2 Sum[sx[S, n, k] , {k, n}] - 

Sum[sz[S, n, k].sz[S, n, Mod [k+1 ,n, 1] ] , {k, n}] 


3.4.2 asymptotic ground states 

The asymptotic ground states for b - 0 and b — * +oo mentioned above are all product 
states of the form \if/) = | Q)® N where \0) is the state of a single spin. We form an N- 
particle tensor product state of such single-spin states with 

in[235] :=productstate [state_?VectorQ , n_Integer/;n>=l] : = 

Flatten [KroneckerProduct @@ Table [state, -fn}] ] 


in accordance with In [160] on page 41. 

The particular single-spin states|f z ), II*), |f z ), |[ z ) we will be using are 

In[236] : = xup [S_?SpinQ] : = 

2~ (-S) Table [Sqrt [Binomial [2S ,M+S] ] ,{M,S, -S, -1}] 
in[237] :=xdn [S_?SpinQ] : = 

2~ (-S) Table [(-1) “ (M+S) Sqrt [Binomial [2S,M+S]] , 
{M,S,-S,-1}] 

in[238] : = zup [S_?SpinQ] := SparseArray [1 -> 1, 2S+1] 
in[239] : = zdn [S_?SpinQ] := SparseArray [-1 -> 1, 2S+1] 
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We can check that these are correct with 


in [ 240 ] :=Table [sx [S] .xup [S] == S*xup[S], {S, 0, 4, 1/2}] 

out [ 240 ] = {True , True, True, True, True, True, True, True, True} 

in [ 241 ] :=Table [sx[S] ,xdn[S] == -S*xdn[S], {S, 0, 4, 1/2}] 

out [241] = {True , True, True, True, True, True, True, True, True} 

in[ 242 ] :=Table [sz [S] . zup [S] == S*zup[S], {S, 0, 4, 1/2}] 

out [ 242 ] = {True , True, True, True, True, True, True, True, True} 

in [ 243 ] :=Table [sz [S] . zdn [S] == -S*zdn[S], {S, 0, 4, 1/2}] 

out [ 243 ] = {True , True, True, True, True, True, True, True, True} 


3.4.3 Hamiltonian diagonalizadon 

We find the m lowest-energy eigenstates of this Hamiltonian with the procedures 
described in subsection 1.7.4: for example, with S = 1/2 and N — 20, 

In [ 244 ] : = With[{S = 1/2, n = 20}, 

(* Hamiltonian *) 
h[b_] = H [S , n, b] ; 

(* two degenerate ground states for b=0 *) 
gsOup = productstate [zup [S] , n] ; 
gsOdn = productstate [zdn [S] , n] ; 

(* ground state for b=+Infinity *) 
gsplusinf = productstate [xup [S] , n] ; 

(* ground state for b=-Infinity *) 
gsminusinf = productstate [xdn[S] , n] ; 

(* numerically calculate lowest m states *) 

Clear [gs] ; 

gs [b_?NumericQ, m_Integer /; m>=l] := gs[b, m] = 
-Eigensystem[-h[N[b]] , m, 

Method -> {"Arnoldi", "Criteria" -> "RealPart", 
Maxlterations -> 10~6}]] 


Comments: 

• gsOup = | i//o[ > and gsOdn = (i//o[) are the exact degenerate ground state wave- 
functions for b — 0; gsplusinf = \y/ +00 ) and gsminusinf = |i^-oo) are the ex- 
act nondegenerate ground state wavefunctions for b = ±oo. 

• The function gs, which calculates the m lowest-lying eigenstates of the Hamil- 
tonian, remembers its calculated values (see subsection 1.6.3): this is impor- 
tant here because such eigenstate calculations can take a long time. 

• The function gs numerically calculates the eigenvalues using h [N [b] ] as a 
Hamiltonian, which ensures that the Hamiltonian contains floating-point machine- 
precision numbers instead of exact numbers in case b is given as an exact 
number. Calculating the eigenvalues and eigenvectors of a matrix of exact 
numbers takes extremely long (please try!). 

• When the ground state is degenerate, which happens here for b = 0, the Arnoldi 
algorithm has some difficulty finding the correct degeneracy. This means that 
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gs [0 , 2] may return two non-degenerate eigenstates instead of the (correct) 
two degenerate ground states (please try for N -15 and N = 20) . This is a well- 
known problem that can be circumvented by calculating more eigenstates: try 
gs [0,10] and check that the two lowest energy eigenvalues are the same. 


• A problem involving N spin-S systems leads to matrices of size (2S + 1 ) ,v x 
(2S + 1 ) v . This scaling quickly becomes very problematic and is at the cen- 
ter of why quantum mechanics is difficult. Imagine a system composed of 
N - 1000 spins S- 1/2: its state vector is a list of 2 1000 = 1.07 x 10 301 complex 
numbers! Comparing this to the fact that there are only about 10 8 ° particles in 
the universe, we conclude that such a state vector (wavefunction) could never 
be written down and therefore the Hilbert space method of quantum mechan- 
ics we are using here is fundamentally flawed. But as this is an introductory 
course, we will stick to this classical matrix-mechanics formalism and let the 
computer bear the weight of its complexity. Keep in mind, though, that this is 
not a viable strategy for large systems, as each doubling of computer capacity 
only allows us to add a single spin to the system, which, using Moore’s law, 
allows us to add one spin every two years. 4 

There are alternative formulations of quantum mechanics, notably the path- 
integral formalism, which partly circumvent this problem; but the computa- 
tional difficulty is not eliminated, it is merely shifted. Modern developments 
such as matrix-product states try to limit the accessible Hilbert space by lim- 
iting calculations to a subspace where the entanglement between particles is 
bounded. This makes sense since almost all states of the huge Hilbert space 
are so complex and carry such complicated quantum- mechanical entangle- 
ment that (i) they would be extremely difficult to generate with realistic Hamil- 
tonians, and (ii) they would decohere within very short time. 


3.4.4 analysis of the ground state 

energy gap 

Much of the behavior of our Ising spin chain can be seen in a plot of the energy 
gap, which is the energy difference between the ground state and the first excited 
state. With m = 2 we calculate the two lowest-lying energy levels and plot their energy 
difference as a function of the parameter b: 


in [ 245 ] :=With[-[bmax = 3, db = 1/64, m = 2 } , 

ListLinePlot [Table [{b, gs [b ,m] [[1,2]] -gs [b,m] [[1,1]]}, 
{b , -bmax , bmax , db}] ] ] 


4 Moore’s law is the observation that over the history of computing hardware, the number of transistors 
on integrated circuits doubles approximately every two years. From http : / /en. Wikipedia, org/wiki/ 
Moore ’ s_law 
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b 


Even in this small 20-spin simulation we can see that this gap is approximately 


Ei~Eq~ 



if \b\ < 1, 
if \b\ > 1. 


(3.31) 


This observation of a qualitative change in the excitation gap suggests that at b — ± 1 
the system undergoes a quantum phase transition (i.e., a phase transition induced 
by quantum fluctuations instead of thermal fluctuations) . We note that the gap of 
Equation (3.31) is independent of the particle number N and is therefore a global 
property of the Ising spin ring, not a property of each individual spin (in which case 
it would scale with N). 


overlap with asymptotic wavefunctions 

Once a ground state wavefunction \y/b) has been calculated, we compute its overlap 
with the asymptotically known wavefunctions with scalar products. Notice that for 
b — 0 we calculate the scalar products with the wavefunctions lyo| ^l^ 01> as they are 
the approximate ground states for \b\ « 1. 


in[246] :=With [{bmax = 3, db = 1/64, m = 2}, 

ListLinePlot [ 

Table [{{b, Abs [gsminusinf . gs [b ,m] [[2,1]]] "2}, 

{b , Abs [gsplusinf . gs [b , m] [ [2 , 1] ] ] ~2> , 

{b, Abs [ ( (gsOup-gsOdn) /Sqrt [2] ) . gs [b,m] [ [2 , 1] ] ] ~2>, 
-Cb, Abs [( (gsOup+gsOdn) /Sqrt [2] ) . gs [b,m] [ [2 , 1] ] ] ~2>, 
{b, Abs [( (gsOup-gsOdn) /Sqrt [2] ) .gs [b,m] [[2, 1]] ] ~2 + 
Abs [( (gsOup+gsOdn) /Sqrt [2] ) .gs [b,m] [[2,1]]] ~2» , 
{b, -bmax, bmax, db}] // Transpose]] 
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b 

Observations: 

• The overlap \(if'b\V / -oo)\ 2 (red) approaches 1 as b — *• -oo. 

• The overlap Ki//j,[i// +00 )| 2 (green) approaches 1 as +oo. 

• The overlap | (y/ b \ | ( C y an j j s mostly negligible. 

• The overlap ^oO^oO | ( oran g e ) approaches 1 as b — ► 0. 

• The sum of these last two, |(y b | | + \{yr b \ | - \{y/ b \y/ 0] )\ 2 + 

\(y/b IV'oj)! 2 (black), approaches 1 as b — ► 0 and is less prone to numerical 
noise. 

• If you redo this calculation with an odd number of spins, you may find differ- 
ent overlaps with the asymptotic wavefunctions. Their sum, how- 

ever, drawn in black, should be insensitive to the parity of N. 

• For \b\ < 0.2 the excitation gap (see above) is so small that the calculated 
ground-state eigenvector is no longer truly the ground state but becomes mixed 
with the first excited state due to numerical inaccuracies. This leads to the 
jumps in the orange and cyan curves (notice, however, that their sum, shown 
in black, is stable). If you redo this calculation with larger values for m, you 
may get better results. 

magnetization 

Studying the ground state directly is of limited use because of the large amount of in- 
formation contained in its numerical representation. We gain more insight by study- 
ing specific observables, for example the magnetization (S®). We add the following 
definition to the With [] clause in In [244] (page 62): 

(* spin components expectation values *) 

Clear [mx , my ,mz] ; 

mx [b_?NumericQ , m_Integer /; m >= 1, k_Integer] := 
mx [b , m, k] = With[{g = gs [b,m] [ [2, 1] ] } , 

Re [Conjugate [g] . (sx[S, n, Mod[k, n, 1] ] . g) ] ] ; 
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my [b_?NumericQ , m_Integer /; m >= 1, k_Integer] := 
my [b , m, k] = With[{g = gs [b,m] [ [2, 1] ]} , 

Re [Conjugate [g] . (sy[S, n, Mod[k, n, l]].g)]]; 
mz [b_?NumericQ , m_Integer /; m >= 1, k_Integer] := 
mz[b, m, k] = With[{g = gs [b ,m] [ [2, 1] ] > , 

Re [Conjugate [g] . (sz[S, n, Mod[k, n, 1] ] . g) ] ] ; 


In our transverse Ising model only the x-component of the magnetization is nonzero. 
Due to the translational symmetry of the system we can look at the magnetization of 
any spin, for example the first one (k - 1): m x (b) (blue) and m z {b) (yellow, non-zero 
due to numerical inaccuracies) 



b 


We see that in the phases of large \b\, the spins are almost entirely polarized, while 
in the phase | b\ < 1 the x-magnetization is roughly proportional to b. 


correlations 

Another very useful observable is the spin-spin correlation function 

*( fc) *(*') *(*) -Ale') 

C k ,k' - (S S )-{S )-(S >. (3.32) 

For S- 1/2 this correlation function has the following known values: 

• — | <C kjk r<+\ 

• Ck,k' = _ § if ih e two spins form a singlet, i.e., if they are in the joint state 

. Remember that the spin monogamy theorem states that if spins k 
and k' form a singlet, then both must be uncorrelated with all other spins in 
the system. 

* Cfc.fc' = 0 for uncorrelated spins. 

* Ck,k' - + j for parallel spins, for example in the Dicke states Iff), , lil), 

or the joint states and . 

We add the following definition to the With [] clause in In [244] (page 62): 
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(* spin-spin correlation operator *) 

Clear [Cop] ; 

Cop [kl_Integer , k2_Integer] := Cop[kl, k2] = 

With [-Cql = Mod[kl ,n, 1] , q2 = Mod[k2,n, 1]}, 

sx [S ,n, ql] . sx [S ,n, q2] + sy [S,n,ql] . sy [S,n,q2] 

+ sz [S ,n,ql] . sz [S,n,q2] ] ; 

(* spin-spin correlations *) 

Clear [c] ; 

c [b_?NumericQ ,m_Integer/ ;m>=l , -[kl_Integer ,k2_Integer}] : = 
c [b,m,{kl ,k2>] = With[{g = gs [b,m] [ [2 , 1] ] > , 

Re [Conjugate [g] . (Cop [kl ,k2] . g) ] - (mx [b ,m,kl] *mx [b ,m,k2] 
+my [b ,m,kl] *my [b ,m,k2] +mz [b ,m,kl] *mz [b ,m,k2] )] ; 


Since our spin ring is translationally invariant, we can simply plot Cg = Cij+g: for 
N — 20 and S — 1 ... 10 (top to bottom), 



b 

Observations: 

• The spins are maximally correlated (C = +|) for b = 0, in the ferromagnetic 
phase. They are all either pointing up or pointing down, so each spin is cor- 
related with each other spin; keep in mind that the magnetization vanishes at 
the same time (page 66). It is only the spin-spin interactions which correlate 
the spins’ directions and therefore their fluctuations. 

• The spins are uncorrelated (C — *• 0) for b — ► ±oo, in the paramagnetic phases. 
They are all pointing in the +x direction for b » 1 or in the -x direction for 
b « -1, but they are doing so in an independent way and would keep pointing 
in that direction even if the spin-spin interactions were switched off. This 
means that the fluctuations of the spins’ directions are uncorrelated. 

entropy of entanglement 

We know now that in the limits b — ■ ±oo the spins are uncorrelated and polarized, 
while close to b - 0 they are maximally correlated but unpolarized. Here we quantify 
these correlations with the entropy of entanglement, which measures the entangle- 
ment of a single spin with the rest of the spin chain. 
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Assume our quantum-mechanical system can be split into two parts, A and B. 
In our case, A is the first spin, and B is the rest of the spin ring; but what follows is 
much more general. Let {| i^>} be a basis set for the description of part A, and {| is)} 
a basis set for the description of part B. In all generality the density matrix of the 
whole system can be expressed as 

pab= £ c i Aij B, i ' A ,j' B \ i A)(i' A \®\iB)(i' B \ ( 3 - 33 ) 

In the case of a pure state \y/), for example when we calculate a ground state of 
our Ising ring, the density matrix is pab = [i/ / )h/ / l- Since we can represent |i//> = 
Li a, j B (puje^A) ® I ]b), *e density matrix is 


Pab - WW = 


£ (pi A ,j B \i A )®\jB) 

£ y 

. i A ,]B 

J'a’J'b 


= £ r }\iA)(i' A \®\jB)(j' B \, (3.34) 


Ia,1' a ,]b,]' b 


which is of the form of Equation (3.33). 

We define the reduced density matrix p A of subsystem A by eliminating subsys- 
tem B through a partial trace : 


p A = Ti B p AB =£0bIPabI7b> = £ c i A jBA'j' liA){i 'A\® 0bI7b>0bIJb> 

J'b i A ,i' A J B ,j' B ,j' B 

= £ Ci A ,iB,i' A ,j B \iA)<i' A l (3.35) 
i A .i' A ,j B 

This density operator only acts on subsystem A and describes its behavior under the 
assumption that we have no access to observables on subsystem B. For a pure state 
[Equation (3.34)], c i Al j n ,i' A ,j' B ~ < t ) ‘Ajn c P*p y < an h the reduced density matrix is 

Pa = £ <!>i A ,jB<t>*i> jB MaXi'aI (3.36) 

i A ,i' A .j B A ’ 


The entropy of entanglement is defined as the von Neumann entropy of the re- 
duced density matrix, 


Sab = -Tr(p A log 2 p A ) = -£A,log 2 A ; (3.37) 

i 

where the A ; are the eigenvalues of p a- Care must be taken with the case A; = 0: we 
find lim^o^-loga^ - 0. 

In Mathematica we define the entanglement entropy of the last spin with the rest 
of the spin ring as follows: 

in [ 247 ] := s [0] = 0; s[x_] = -x Log [2 , x] ; 

in [248] :=EE [S_?SpinQ , psi_] := Module [{g, rhoA}, 

(* group the wavefunction coefficients into lists, *) 

(* such that the state of the last spin is the same *) 

(* for all basis states within that list *) 

g = Transpose [Partition [psi , 2S+1]]; 
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(* compute the single-spin reduced density matrix *) 

(* of the last spin *) 

rhoA = Conjugate [g] . Transpose [g] ; 

(* compute entropy of entanglement *) 

Total [s /<§ Re [Eigenvalues [rhoA] ]] ] 


Observations: 


• Entanglement entropies of the known asymptotic ground states: 


1 

in [249] : =EE [1/2 , (gsOup+gsOdn) /Sqrt [2] ] 

2 

Out [249] = 1 

3 

in [250] : =EE [1/2 , (gsOup-gsOdn) /Sqrt [2] ] 

4 

Out [250] = 1 

5 

in[25i] :=EE [1/2 , gsplusinf] 

6 

Out [251] = 0 

7 

in[252] :=EE [1/2 , gsminusinf] 

8 

Out [252] = 0 


• Entanglement entropy as a function of b: again the calculation is numerically 
difficult around b =; 0 because of the quasi-degeneracy. 


2 

3 


in[253] :=With [{bmax = 3, db = 1/64, m = 2} , 

ListLinePlot [Table [{b , EE[l/2, gs [b ,m] [ [2 , 1] ] ] } , 
{b, -bmax, bmax, db}] , PlotRange -> {0, 1}]] 



b 

Notice that the quantum phase transition is not visible in this plot. 

3.4.5 exercises 

Q3.14 Show that the single-spin states of Equation (3.28), implemented in In [236] 
and In [237] , are indeed eigenstates of S x with eigenvalues +S. 

Q3.15 For S = 1/2, what is the largest value of N for which you can calculate the 
ground-state wavefunction of the transverse Ising model at the critical point 
b= 1? 
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Q3.16 Study the transverse Ising model with S = 1: 

1 . At which values of b do you find quantum phase transitions? 

2. Characterize the ground state in terms of magnetization, spin-spin cor- 


1 . Guess the shape of the ground-state wavefunctions for b± oo [notice that 
the first term in the Hamiltonian of Equation (3.38) is in the z-direction!] 
and compare to the numerical calculations. 

2. At which values of b do you find quantum phase transitions? 

3. Characterize the ground state in terms of magnetization, spin-spin cor- 
relations, and entanglement entropy. 


Q3.18 Do the same calculation for S — 1/2 with the Heisenberg-interaction Hamilto- 
nian 


1 . Guess the shape of the ground-state wavefunctions for b± oo [notice that 
the first term in the Hamiltonian of Equation (3.38) is in the z-direction!] 
and compare to the numerical calculations. 

2. What is the ground-state degeneracy for b- 0? 

3. At which values of b do you find quantum phase transitions? 

4. Characterize the ground state in terms of magnetization, spin-spin cor- 
relations, and entanglement entropy. 


Q3.19 Consider two spin- 1/2 particles in the triplet state | if/) = ] ] ). SubsystemAis 
the first spin, and subsystem B is the second spin. 

1. What is the density matrix / oab of this system? 

2. What is the reduced density matrix pa of subsystem A (the first spin)? Is 
this a pure state? If yes, what state? 

3. What is the reduced density matrix ps of subsystem B (the second spin)? 
Is this a pure state? If yes, what state? 

4. Calculate the von Neumann entropies of pab , p a , and p / j . 

Q3.20 Consider two spin-1/2 particles in the singlet state \y/) - ' ' • Subsystem 

A is the first spin, and subsystem B is the second spin. 

1. What is the density matrix p ab of this system? 

2. What is the reduced density matrix pa of subsystem A (the first spin)? Is 
this a pure state? If yes, what state? 

3. What is the reduced density matrix ps of subsystem B (the second spin)? 
Is this a pure state? If yes, what state? 

4. Calculate the von Neumann entropies of p a b , P a . and ps- 


relations, and entanglement entropy. 
Q3.17 Study the transverse XY model for S = 1/2: 



(3.38) 





(3.39) 


Chapter 4 

real-space systems 


4. 1 one particle in one dimension 


One-dimensional single-particle systems are governed by Hamiltonians of the form 




h 2 d 2 

~2ine^* nxl 


(4.1) 


The system’s behavior is determined by the mass m and the potential V(x). 

In what follows we restrict the freedom of the particle to a domain x e Q = [0 ,a], 
where a can be very large in order to approximately describe quasi-infinite systems. 
This assumes the potential to be 


{ oo for x < 0 

W(x) forO <x<a (4.2) 

oo for x > a 

This restriction is necessary in order to achieve a finite representation of the system 
in a computer. 

4.1.1 basis functions 

The Hilbert space of this particle consists of all square-integrable (L 2 ) and differ- 
entiable wavefunctions with support in Q. For each ket \ij/) in this Hilbert space 
we define the wavefunction y/(x) - (x\y/) in terms of the “position basis” {\x)\ x€ (i, 
which satisfies the completeness relation / 0 ” |x)(x|dx = In- 1 The scalar product in Q 
is defined as 

r pa pa pa 

<V / ll> = <Vl / U){x\&x lx)- (y/\x)(x\x)dx= / iy*(x)xWdx. (4.3) 
[JO Jo JO 

As usual we need a set of basis functions {| i)},* to describe this system. There 
are many possible choices of basis functions. The position basis {|x)} xe n is ill suited 

J To be exact, the position basis set {|x)} xe Q spans a space that is much larger than the Hilbert space 
of square-integrable smooth functions used in quantum mechanics. This can be seen by noting that this 
basis set has an uncountable number of elements, while the dimension of the Hilbert space in question 
is only countably infinite [see Equation (4.4) for a countably infinite basis set] . For example, the state 
Xxe(nnQ) I*) i s not a valid quantum-mechanical state (it is too pathological), yet it can be expressed in 
this position basis. 
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for this task, since its elements are singular and therefore difficult to represent in a 
computer. The most generally useful ones are the momentum basis and the finite- 
resolution position basis, which we will look at in turn, and which will be shown to 
be related to each other by a type-I discrete sine transform. 


momentum basis 

The simplest one-dimensional quantum-mechanical system of the type of Equa- 
tion (4.1) is the infinite square well with W(x) = 0. Its energy eigenstates are 


(x\n) =<pn(x) = 



(4.4) 


for n = 1,2,3, ..., with eigen-energies 


_ n 2 n 2 k 2 
n 2 ma 2 


(4.5) 


We know from the Sturm-Liouville theorem that these functions form a complete 
set (see Q2.2 on page 35); further, we can use Mathematica to show that they are 
ortho - normalized: 


1 in [ 254 ] :=phi [a_ , n_ , x_] = Sqrt [2/a] Sin[n Pi x/a] ; 

2 in [ 255 ] :=Table [Integrate [phi [a, nl ,x] *phi [a, n2,x] , {x, 0, a}], 

3 {nl, 0, 10}, {n2, 0, 10}] // MatrixForm 


( \2 2 


p 2 \n) = 


n 2 7t 2 H 2 

2 !”>■ 

a A 


(4.6) 


This makes the kinetic-energy operator = p 2 / [2m) diagonal in this basis: (n\J{Y m \n') 
E n Snn'- However, in general the potential energy, and most other important terms 
which will appear later, are difficult to express in this momentum basis. 

The momentum basis of Equation (4.4) contains a countably infinite number of 
basis functions. In practical calculations, we restrict the computational basis to ne 
{1 . . . n m axl> which means that we only consider physical phenomena with excitation 
energies below £„ max = 


finite-resolution position basis 

Given an energy-limited momentum basis set {| we define a set of n mea equally- 
spaced points 

x; = ax - (4.7) 

ttmax + 1 

for j e {1 . . . ttmaxl- We then define a new basis set as the closest possible representa- 
tions of delta-functions at these points: 


l}> = 


/ 


a 

^max + 1 


^max 

£ (f>n[Xj)\n). 


72=1 


(4.8) 
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The spatial wavefunctions of these basis kets are 


(x\j) = R j (x) = J Y, (t>n(Xj)(p n {x). (4.9) 

V ^max + 1 n-l 

Here is an example of what these position-basis functions look like for n max =10: 





This new basis set is also ortho-normal, (j\j f ) = 6 jjr, and it is strongly local in the 
sense that only the basis function i) j (x) is nonzero at xj, while all others vanish: 

(Xj'lj) = 9j (xy) = 5jj' X \/ n ' T ' a ^ + 1 ■ (4.10) 

We define these basis functions in Mathematica with 

in [256] := nmax = 10; 

in [ 257 ] :=xx [a_ , j _] = a j/(nmax+l); 

in[258] :=theta[a_, j_, x_] = 

Sqrt [a/(nmax+l)l Sumfphi [a,n,xx [a, j] ] phi[a,n,x], 

{n, 1 , nmax}]; 


Since the basis function 9 j{x) is the only one which is nonzero at Xj, and it is 
close to zero everywhere else (exactly zero at the xy^j), we can usually make two 
approximations: 

• If a wavefunction is given in the position basis, \y/) = Vj\j), then by 

Equation (4.10) the wavefunction is known at the grid points, if/[xj) = vj x 

^/ «max+i _ |p|j s a ]] ows for very easy plotting of wavefunctions and densities by 
linearly interpolating between these grid points: 

1 In[259] : = ListLinePlot [ 

2 Transpose [{Table [xx [j] , {j , 1, nmax}], 

3 (nmax+l)/a * Abs[v]~2}]] 


By the truncation of the basis at n max , the wavefunction has no frequency 
components faster than one half-wave per grid-point spacing, and therefore 
we can be sure that this linear interpolation is a reasonably accurate represen- 
tation of the full density | (x|i//> | 2 , in particular as « rnax — ► oo. 
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• If a potential energy function W (x) varies smoothly over length scales of jty+i - 
Xj - "/("max + 1), then the matrix elements of this potential energy in the po- 
sition basis are approximately diagonal, 


% = (j\V\f) = (j\ [ jjx><x|dx V jT a |jc / ><x'|d*' 


I /> 


= J d xj dx' (j\x)(x\V\x'){x'\j') = J dxdj(x)W(x)'Oji(x ) 

J fCl 

dxd*(x)i h'{x) = 8 jjiW(Xj), (4.11) 
0 J 


where we have used Equation (4.2) and the fact that {x\ V\x') = W(x)8{x- x') 
since the potential is diagonal in the position basis, as well as the approximate 
locality of Rj[x) around xj implying W{x)i)j(x) ~ W(Xj)f)j{x). This is a mas- 
sive simplification compared to the explicit evaluation of potential integrals 
for each specific potential energy function. 


conversion between basis sets 

Within the approximation of a truncation at maximum energy E „ m;ly , we can express 
any wavefunction |i//> in both basis sets of Equation (4.4) and Equation (4.9): 

ttmax TZmax 

W) = Y u n | n) - Y Vj\j) (4-12) 

n = 1 7=1 

Inserting the definition of Equation (4.8) into Equation (4.12) we find 


^max ^max 

Y u n\n)= Y v i 

n= 1 7=1 


(2 '‘max 

— Y <Pn’(Xj)\ri) 

W-max + 1 n'-l 


= Y 

n'-l 


— Y Vj<Pn'(Xj) 

n max 1 j- \ 


I n!) (4.13) 


and therefore, since the basis set {| n) j is ortho-normalized, 


u n — 


/ 


a 

^max + 1 


^max ^max 

Y v .i ( Pn(Xj) = Y X njVj 
7=1 7=1 


with the basis conversion coefficients 


(4.14) 


nnj \ 

"max + 1 i 

(4.15) 

The inverse transformation is found from \n) - (j\n)\j) inserted into Equa- 

tion (4.12), giving 

^max 

Vj— ^ X n jU n (4.16) 

n= 1 

in terms of the same coefficients of Equation (4.15). Thus the transformations relat- 
ing the vectors u (with components u n ) and v (with components vj) are v = X ■ u 
and u-Xv in terms of the same symmetric matrix X with coefficients X n j. 

We could calculate these coefficients with 


X n j - 


c + 1 


<M*/) : 


: + 1 


2 

— sin 
a 


( nnxi \ 

in hH : 


r + 1 


- sm 
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in [ 260 ] :=X = Table [Sqrt [2/ (nmax+1 ) ] Sin[Pi*n*j/(nmax+l)] , 
{n, 1, nmax}, {j , 1, nmax}-] // N; 


but this is not very efficient, especially for large n max . 

It turns out that Equation (4.14) and Equation (4.16) relate the vectors u and 
v by a type-I discrete sine transform (DST-I), which Mathematica can evaluate very 
efficiently via a fast Fourier transform/ Since the DST-I is its own inverse, we can 
use 

in[26i] :=v = Four ierDST [u , 1]; 
in [262] :=u = Four ierDST [v , 1]; 


to effect such conversions. We will see a very useful application of this transforma- 
tion when we study the time-dependent behavior of a particle in a potential (“split- 
step method”, subsection 4.1.3). 

The matrix X is calculated most efficiently by repeated calls to the DST-I func- 
tion: 

in [263] : = SparseIdentityMatrix [n_] : = 

SparseArray [Band[{l , 1}] ->1 , {n,n}-] 
in [264] :=X = Four ierDST [# , 1] & /<§ SparseldentityMatrix [nmax] ; 


This matrix notation of the basis transformation is useful for converting operator 
representations between the basis sets: the momentum representation U and the 
position representation V of the same operator satisfy 

In [265] : = V = X.U.X; 

In[266] :=U = X.V.X; 


In practice, as above we can convert operators between the position and momentum 
representation with a two-dimensional type-I discrete sine transform: 

in [267] :=V = Four ierDST [U , 1]; 
in[268] :=U = Four ierDST [V , 1]; 


This easy conversion is very useful for the construction of the matrix representations 
of Hamiltonian operators, since the kinetic energy is diagonal in the momentum 
basis, Equation (4.6), while the potential energy operator is approximately diagonal 
in the position basis, Equation (4.11). 


special case: the kinetic energy operator 


The representation of the kinetic energy operator can be calculated exactly with the 
description given above. Using Equation (4.6), the kinetic Hamiltonian is 


p z 1 

3P, . - r ~ 

•sikin — _ ~ „ 

2m 2m 


n max 
n= 1 


£ I n'){ri | 

n'= 1 


^max 

= Ei £ n z \n){n\, (4.17) 

n-l 


2 fc2 

where E i = [see Equation (4.5)]. In Mathematica, we define the kinetic energy 

operator in the momentum basis as 

2 see http://en.wikipedia.org/wiki/Fast_Fourier_transform 
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in [269] :=HkinM = H ~ 2 Pi~2/(2 m a~2) * 

SparseArray [Band [{1,1}-] ->Range [nmax] "2] 


From this we can calculate the representation in the finite-resolution position basis 
with 

in [ 270 ] :=HkinP = FourierDST[HkinM, 1]; 


However, for large n max it is often acceptable to use the following approximation: 


O’l^kinl/) * £l x 


n max( n max+2) j |- j _ jl 

(—1)7—/ x 2 wy if y ^ /. 


(4.18) 


We will not be using this approximation in what follows, as the basis-set conversion 
through type-I discrete sine transforms is usually sufficiently efficient. 


4.1.2 example: square well with bottom step 

A simple example you may remember from quantum-mechanics class is a particle 
moving in the one-dimensional potential given by Equation (4.2) with 


where we assume 
fashion. 


W(x) = 


Wo 

0 


if x< | 
if x> f 


(4.19) 


Wo > 0 for simplicity; the case Wq < 0 is solved in the exact same 


> 



x / a 


analytic solution for 0 < E < W () 

The potential of Equation (4.19) is so simple that we can find the eigenstates of this 
particle analytically. Since the potential is piecewise flat, we know that the energy 
eigenstates must be (hyperbolic) sine and cosine functions with piecewise constant 
wavelengths. In order to find these wavelengths we set 

r x 

W\ (x) = A sinh \k\Ti — 
y a 


a 

for 0 < x < - 
2 
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which satisfy i// j (0) = i// 2 (a) = 0 to match the boundary conditions where the poten- 
tial becomes infinite. We assume that k\,k 2 > 0. 

The conditions for matching these two pieces of the wavefunction are i/r 1 ( § ) = 
1^2 ( § ) and i/Zj ( | ) = y/' 2 ( | ) , from which we find the condition 


nk\ nk 2 

k\ coth -~ko cot . 

2 2 


(4.21) 


The time-independent Schrodinger equation further supplies the energy condition 


E= W 0 


h 1 2 3 n 2 k\ 
2 mci 2 


h 2 n 2 k 2 
2 ma 2 


(4.22) 


Since we have assumed that 0 < E < Wo we find from this that k\- \J El - k 2 in terms 
of the dimensionless parameter 


Wq 2ma 2 Wo 
Ei n 2 H 2 


(4.23) 


We notice that the entire problem only depends on this one dimensionless param- 
eter Q, and not on the individual values of m, a, and Wq: the effort of making the 
problem dimensionless has paid off by significantly reducing the number of param- 
eters that we need to study. The resulting eigenvalue equation 


Q. - k 2 coth — =—k 2 cot — — . (4.24) 

thus depends only on one parameter Q, and can be solved graphically for k 2 in the 
range 0 < k 2 < Vo.. 

For O < 1.66809 there is no solution for k 2 , meaning that the ground state has 
energy E > Wq. As a numerical example, for El - 2 we plot the left-hand side of 
Equation (4.24) in blue and the right-hand side in red: 

1 in[27i] :=With [{Omega = 2}, 

2 Plot [{Sqrt [0mega-k2~2] Coth [Pi Sqrt [0mega-k2~2] /2] , 

3 -k2 Cot [Pi k2/2]}, {k2, 0, Sqrt [Omega] }] ] 



0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 


k 2 
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We find a single solution for the ground state at k 2 = 1.32884 numerically with 

in [ 272 ] := s = With [{Omega = 2}, 

FindRoot [Sqrt [0mega-k2~2] Coth [Pi Sqrt [0mega-k2~2] /2] 
== -k2 Cot [Pi k2/2] , {k2, 1}]] 

Out [ 272 ] = {k2 -> 1.32884> 


Notice that the result is given as a list of replacement rules (with the -> operator). 
You can extract the value of k 2 with 

In [273] : = k2 /. S 
out [273] = 1 . 32884 


and calculate the value of Aq with 

in[274] :=With [{Omega = 2}, 

in [ 275 ] := Sqrt [0mega-k2~2] /. s] 

out [275] = 0.48392 


We can plot the result with (assuming a — 1) 

in [276] :=With [{kl=0 . 4839202839634602 , k2=l . 3288420368007343 , 

A=1 . 6088142613650431 , B=1 . 5458263302568298} , 
psi0[x_] = Piecewise [{{A Sinh[kl Pi x] , 0<=x<=l/2}, 

{B Sin [k2 Pi (1-x)] , 1/2<x<=1}>] ; 
Plot [psiO [x] , {x, 0, 1}, Exclusions->None] ] 



We will be using this wavefunction psiO [x] below for a comparison with numerical 
calculations. 

For E > Wq the same calculation must be re-done with i//| (x) = i4sin(Aq xla). 
The algebra is very similar, and the results do not teach us anything further for this 
course. 

exercises 

Q4. 1 Find and plot the ground state for Q = 1000. What is the probability to find the 
particle in the left half of the potential well? 
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numerical solution (I): momentum basis 


We first search for the ground state of the step-well in the momentum basis. The 
matrix elements of the kinetic energy are diagonal, 


< ll | H ) — — o - x filin’- 

The matrix elements of the potential energy of Equation (4.19) are 


(4.25) 


(n\V\n‘ 


pa 

') = 4>n 

Jo 


(x)W(x)(p n i (x)dx 


2 Wp 

a Jo 


a 

/„' sln ( 


nnx\ 
a I 


sin 


) f 2 sir 
Jo 


= Wo x < 


dx = 2Wol s\n{rmy)sin{n' ny)Ay 
if n - n' 

if n + n' odd (4.26) 
otherwise 


n+n'+l 


n-n ' +1 


This allows us to express the matrix elements of the Hamiltonian H nn > = {n\J€\n!) 

2 fc2 

in units of the energy E\ = : 


a nn‘ 

“eT 


= n8 ntl i + D. x < 


n+n ' + 1 

c-l) z 

n+n' 


n-n ! + 1 
n-n' 


if n = n' 
if n + n' odd 
otherwise 


(4.27) 


with Q = Wq I Ei the same dimensionless parameter as above. In Mathematica, 

in[277] :=h [Qmega_ , n_, n_] = n~2 + Qmega/2; 

in[278] :=h [Omega_ , n_ , np_] /; OddQ [n+np] = Qmega/Pi * 

((-l)~((n+np+l)/2)/ (n+np) - (-1) “ ( (n-np+1) /2)/ (n-np) ) ; 
in[279] :=h [Omega_ , n_ , np_] /; EvenQ [n+np] = 0; 


For a given « max we can now find the Hamiltonian matrix with 
in [280] :=nmax = 10; 

in[28i] :=H [0mega_] = Table [h [Omega, n,np] , {n,l,nmax}, {np , 1 ,nmax}] ; 


and the ground state coefficients with 

in [282] : = Clear [gs] ; 

in[283] :=gs [Omega_?NumericQ] := gs [Omega] = 

-Eigensystem[-H [N [Omega] ] , 1, 

Method -> {"Arnoldi", "Criteria" -> "RealPart", 
Maxlterations -> 10~6}] 


The ground state wavefunction <jc| y„ max ) is then 

in[284] :=psi [Omega_?NumericQ , a_, x_] : = 

gs [Omega] [ [2, 1] ] . Table [phi [a, n, x] , {n, 1, nmax}] 
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For Q = 2 we can calculate the overlap of this numerical ground state with the exact 
one given in In [276] , <Vol7n max >: 

in [285] :=NIntegrate [psiO [x] *psi [2 , 1 ,x] , {x , 0 , 1}] 


This overlap quickly approaches unity as n max increases: 



t^max 

The red line is a least-squares fit to the even values of n max , giving a convergence of 
1 - (Volrnmax) 2 ~ 0.0099 ra” 4 ^ 1 . 

numerical solution (II): mixed basis 

The main difficulty of the first numerical solution above was the evaluation of the 
potential matrix elements of Equation (4.26). For such a simple step potential as 
used here, we were able to find an analytic expression for {n\V\n')\ but for more 
complicated potentials this will not be possible. But we have seen in Equation (4.11) 
that the potential is approximately diagonal in the finite-resolution position basis, 
and we can therefore find an approximate expression for the Hamiltonian matrix 
with a procedure that is independent of the shape of W (x). 

We first calculate the matrix elements of the kinetic energy operator in the mo- 
mentum basis, again in units of E \ : 

in [ 286 ] :=nmax = 10; 

in [287] :=HkinM = SparseArray [Band [{1 , l}] ->Range [nmax] ~2] ; 


Next we convert this operator into the finite-resolution position basis: 
in [ 288 ] :=HkinP = FourierDST [HkinM , 1]; 


The potential energy operator W is expressed approximately in the finite-resolution 
position basis: setting a - 1 for simplicity and using Q = Wq! E\, 

in[ 289 ] :=W [0mega_ , x_] = Piecewise [{{Omega, x<l/2}, 

{Omega/2, x==l/2>, {0, x>l/2»] ; 
in[ 290 ] :=xval = Table [j / (nmax+1) , {j , 1, nmax}]; 
in[ 29 i] :=Wval [0mega_] = W [Omega, #]& /@ xval; 

in[ 292 ] :=HpotP [0mega_] = SparseArray [Band [{1 , 1}] -> Wval [Omega] ] ; 
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The full Hamiltonian in the finite-resolution position basis is 


in[ 293 ] :=HP [Qmega_] = HkinP + HpotP [Omega] ; 


We find the ground state with 


in[ 294 ] : = Clear [gsP] ; 

in[ 295 ] :=gsP [Omega_?NumericQ] := gsP [Omega] = 

-Eigensystem [-HP [N [Omega] ] , 1, 

Method -> {"Arnoldi", "Criteria" -> "RealPart", 
Maxlterations -> 10~6}] 


and, as shown before, this can be plotted simply with 


in [ 296 ] :=With[{0mega=2}, 
gammaP = Join[ 

{{ 0 , 0 }}, 

Transpose [{xval , Sqrt [nmax+1] *gsP [Omega] [[2,1]]}] , 
f{l,0}>] ; 

ListLinePlot [gammaP] ] 


where we have “manually” added the known values y(0) = y (1) = 0 to the list of nu- 
merically calculated wave-function values. 



You can see that even with rc max = 10 grid points this ground-state wavefunction 
(thick blue line) looks remarkably close to the exact one (thin red line, see page 78). 

The wavefunction is calculated by converting to the momentum representation 
as in In [262] and multiplying with the basis functions as in In [284] : 


in[ 297 ] :=psiP [Omega_?NumericQ , a_, x_] : = 

FourierDST [gsP [Omega] [ [2 , 1] ] , 1] . 
Table [phi [a , n , x] , {n , 1 , nmax}] 


As for In [285] the overlap of this numerical wavefunction with the exact one ap- 
proaches unity as n max increases: 
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^max 


The red line is a least-squares fit to the even values of n max , giving a convergence 
of 1 - (Vo\Tn max ) 2 ~ 0.03061 78 w^’ 4889 . This convergence is slower than for the 
momentum-basis calculation since there was an additional approximation involved 
in Equation (4.11). 


exercises 


Q4.2 Find and plot the ground state for Q = 1000, using the approximate numeri- 
cal method. What is the probability to find the particle in the left half of the 
potential well? 


Q4.3 


Calculate the energy levels and energy eigenstates of a particle in a well with 
bottom potential 


W(x) = 




(4.28) 


Compare them to the analytically known eigen-energies and eigenstates of a 
harmonic oscillator. 


Q4.4 With a = 1, take a “noisy” potential W(x) = a n (p n (x) with a n random: 

(a n ) - 0 and (a 2 n ) - n~ 2 . Plot the ground-state density |y(x)| 2 using « max » h, 
for different values of Q. 


4.1.3 dynamics 


Assume again a single particle of mass m moving in a one-dimensional potential, 
with Hamiltonian 




h 2 d 2 
2 


(4.29) 


'^dcirl 


-''f poi 


The motion is again restricted to x e [0, a]. We want to study the time-dependent 
wavefunction y/{x, t) = (x|i//(t)) given in Equation (2.33) on page 39, 


I y{t)) = exp 


i [t- t 0 ) 

h 


Je 


Wtto)). 


(4.30) 


The simplest way of computing this propagation is to express the wavefunction 
and the Hamiltonian in a particular basis and use a matrix exponentiation to find the 
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time dependence of the expansion coefficients of the wavefunction. For example, if 
we use the finite-resolution position basis, we have seen on page 80 how to find the 
matrix representation of the Fiamiltonian, HP. For a given initial wavefunction psiO 
we can then define 

i in[298] :=psi [t_?NumericQ] := MatrixExp [-I*t*HP] . psiO 


where we have changed the units such that the time t = (f- fo) lh is in units of inverse 
energy. If you try this out, you will see that calculating \y/(t)) in this way is not very 
efficient, because the matrix exponentiation is a numerically difficult operation. 

A much more efficient method can be found by first splitting up the Hamiltonian 
as — ^ki n + .#p OI as in Equation (4.29), and then using the Trotter expansion 


e MX+Y) = e §X e AY e §X xe § T lX,[X,Yll + j I [Y,lX.Yll ><e - 


:e iX e AY e iX 


(4.31) 


where the approximation is valid for small A since the neglected terms are of third 
and higher orders in A (notice that there is no second-order term in A!). Setting 
A = - ' U M jf for some large integer M, as well as X = .7Aj )0 i and Y — we find 


M- 


i M 


|i//(f)>= lim \e xje \ \iy(to)) — lim L^kin+^pot) \y/(t 0 )) 


M- 


M—> oo 
A 


1 M 


M 


= lim e 2^pot e 4^i dne2 ^ 0 pot \y/(t 0 )) 


= lim e 2 ^pote A ^un e ^pot e A^ tine A^’p Ot ... e -l^i d n e 2 -^p Ot | 1 «(f 0 )). (4.32) 

M — oo ' ' 

(M - 1) repetitions of e A ^ in 


This can be evaluated very efficiently. We express the potential Hamiltonian in the 
finite-resolution position basis, the kinetic Hamiltonian in the momentum basis, 
and the time -dependent wavefunction in both bases of Equation (4.12): 


n max n max 

= u n (t)\n) = E v W 

n = 1 7=1 

n max 

£ W{xj)\j)(j\ 

7=1 

n max 

•A^kin = E n 2 \n)(n\ 

n = 1 


(4.33) a 

(4.33) b 

(4.33) c 


where we have already expressed all energies as multiples of the square-well ground- 

2 fc2 

state energy E\ = * . The expansion coefficients of the wavefunction are related 

by a discrete Fourier transform, Equation (4.16). 

The matrix exponential of a diagonal matrix is easily found from the Taylor ex- 
pansion: 




ir _ jpk - V __ 
A fc! P ot A lr\ 


k = 0 


k = 0 


k\ 


n max 

k S A fc 

^max 

£ w{ Xj )\j){j\ 

. 7=1 

k = 0 K ■ 

£ wHxj)\mj\ 

. 7=1 


= E 

7=1 


E u wk W 


k = 0 


\j)(j \ = £ e Amx i ] \ (4.34) 
7=1 
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where we have used the integer matrix powers 


£ W[x h )\j 2 ){h\ 

72=1 


£ w(x h )\j k )Uk\ 

ik = 1 


£ W[ Xj )\j){j\ 
7=i 


£ w{x h )\h){h\ 

L7i=i 

^max w max ^max 

= £ £ ••• £ W(x h )W{x h )---W{Xj k )\j 1 )(j 1 \j 2 )(j2\j3)---(jk-i\jk)(jk\ 

71=1 j 2 =l j k =l 

^max w max ^max 

-££-£ \V[x h )W(x h )..AY(x h )j^j^P.jr-^ h : ,/ s (jk\ 

71=1 7' 2 =1 7*=1 


= £ W*(xj)|j><./|. (4.35) 

7=1 


The action of the potential Hamiltonian thus becomes straightforward: 




£ 

7=1 




l7")(7l 


^max 

£ Vj'{t)\j') 

7'=i 


^max r -1 

= £ |y>, (4.36) 

7=i J 


which is a simple element-by-element multiplication of the coefficients of the wave- 
function with the exponentials of the potential - no matrix operations are required. 
The expansion coefficients (position basis) after propagation with the potential Hamil- 
tonian are therefore 


v'j - 


V P 


(4.37) 


The action of the kinetic Hamiltonian in the momentum representation is found 
in the exactly same way: 




I Vd?)> : 


£ e \n)(n\ 

n=l 



n max r 2 1 

= £ \e x,r u n {t) \\n). 

n= 1 L 


(4.38) 


The expansion coefficients (momentum basis) after propagation with the kinetic 
Hamiltonian are therefore 

u' n = e An2 u n . (4.39) 

We know that a type-I discrete sine transform brings the wavefunction from the 
finite-resolution position basis to the momentum basis of Equation (4.16) . The prop- 
agation under the kinetic Hamiltonian thus consists of 

1. a type-I discrete sine transform to calculate the coefficients uj <-*■ u n , 

2. an element-by-element multiplication, Equation (4.39), to find the coefficients 
U n >— u' n , 

3. and a second type-I discrete sine transform to calculate the coefficients u' n =-» 
v 'p 

Here we assemble all these pieces into a program which propagates a state \y/(to)) 
given as a coefficient vector V in the finite-resolution basis forward in time to t ~ 
to + A t with M time steps. We assume a — 1 and E\ - 1. 
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in[299] :=HpotP = SparseArray [Band [{1 , 1}] ->Wval] ; 
in [ 300 ] :=HkinM = SparseArray [Band [-C 1 , 1}] ->Range [nmax] ~2] ; 
in[30i] :=HkinP = FourierDST[HkinM, 1]; 
in [302] :=HP = HkinP + HpotP; 

in[303] :=propExact [Dt_?NumericQ , psi_?(VectorQ[#,NumericQ]&)] : = 
MatrixExp [-I*Dt*HP] .psi 

in [304] :=propApprox [Dt_?NumericQ , M_Integer/ ;M>=2, 

psi_?(VectorQ [#,NumericQ)]&] := 

Module [{Ke ,Pe2 ,Pe ,psil , psi2 , propKin, propPot} , 

(* compute exponentials *) 

Ke = Exp [-I*Dt/M*Table [n~2 , {n, 1 , nmax}] ] ; 

Pe2 = Exp [-I/2*Dt/M*Wval] ; 

Pe = Pe2~2; 

(* propagate with the potential operator *) 
propPot [p_] := Pe * p; 

(* propagate with the kinetic operator *) 
propKin[p_] := FourierDST [Ke*FourierDST [p , 1] , 1] ; 

(* propagation *) 
psil = propKin[Pe2*psi] ; 

psi2 = Nest [propKin [propPot [#]]&, psil, M-l] ; 
Pe2*psi2] 


Notice that there are no basis functions, integrals, etc. involved in this calculation; 
everything is done in terms of the values of the wavefunction on the grid xi . . . x Umax . 
This efficient method is called split-step propagation. 

The Nest command “nests” a function call: for example, Nest [f ,x,3] calcu- 
lates /(/(/(x))))). We use this on line 20 above to repeatedly propagate by the poten- 
tial and kinetic operators. This propagation algorithm can be adapted to calculate 
the wavefunction at all the intermediate times t — to + jjj ( t - to) for m - 1, 2, 3, . . . , M, 
which allows us to follow the evolution of the wavefunction during its time evolu- 
tion. To achieve this we simply replace the Nest command with NestList, which is 
similar to Nest but returns all intermediate results: for example, NestList [f , x , 3] 
calculates the list {x, /(x), /(/(x)) , /(/(/(x)))}. We replace the code above from line 20 
with 

20 
21 
22 


psi2 = NestList [propKin [propPot [#] ] &, psil, M-l]; 
Transpose [{Range [0 , M] *Dt/M , 

Prepend [(Pe2*#) & /<§ psi2, psi]}]] 


exercises 

Q4.5 Convince yourself that the Trotter expansion of Equation (4.3 1) is really neces- 
sary, i.e., that e x+Y e x e Y if X and Y do not commute. Hint: use two concrete 
non-commuting objects X and Y, for example two random 2x2 matrices as 
generated with RandomReal [{0 , 1} , {2 , 2}] . 

Q4.6 Given a particle moving in the range x e [0, 1] with the scaled Hamiltonian 

1 d 2 

— x- — - + Osin(107rx), 

n l dx z 


(4.40) 
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(x-112) 2 

compute its time-dependent wavefunction starting from y/{t - 0) ex e x, 2 e lKX 
with (T = 0.05 and k = 100. Compute (x)(f) for t - 0.. .0.2 using first Q = 0 and 
then O = 1000. 


4.2 Many particles in one dimension: dynamics with the 
non-linear Schrodinger equation 


The advantage of the split-step evolution of Equation (4.32) becomes particularly 
clear if the particle evolves according to the non-linear Hamiltonian 


hr d 2 2 


-^kin 


Xx 


pot 


(4.41) 


Such Hamiltonians can describe the mean-field interactions between N particles 
which are all in wavefunction t//(x), and which are therefore in a joint product wave- 
function y/[x)® N . One particle’s wavefunction \j/{x) (normalized to / [i//(x, f)| 2 dx = 1) 
sees a potential generated by the average density (AT- l)|i/r(x)| 2 of other particles 
with the same wavefunction, usually through collisional interactions. The associ- 
ated non-linear Schrodinger equation is called the Gross-Pitaevskii equation and 
describes the dynamics of Bose-Einstein condensates: 


i h 


dy/{x, t) 
dt 


h^J?_ 

2 m dx 2 


+ V(x) + g\y/{x, f)| 2 


y/(x, t ). 


(4.42) 


The coefficient g - (N - 1) x Airf ^ n as approximates the mean-field 5-wave scattering 
between a particle and the {N- 1) other particles, with 5-wave scattering length a s . 

For any g f 0 there is no solution of the form of Equation (4.30). But the split- 
step method of Equation (4.32) can still be used because the potential is still diag- 
onal in the position representation. We extend the Mathematica code of the pre- 
vious section by modifying the propApprox method to include a non-linear term 
with prefactor g, and do not forget that the wavefunction at grid point xj is yr(X j) = 

in [306] :=propApprox [Dt_?NumericQ , M_Integer/ ;M>=2 , g_?NumericQ, 
psi_?(VectorQ[#,NumericQ]&)] := 

Module [-[Ke,psil ,psi2,propKin,propPot} , 

(* compute exponentials *) 

Ke = Exp [-I*Dt/M*Table [n~2 ,-(n, 1 ,nmax}] ] ; 

(* propagate with the potential operator *) 
propPot [dt_ ,p_] := 

Exp [-I*dt* (Wval+g* (nmax+1) *Abs [p] ~2) ] * p; 

(* propagate with the kinetic operator *) 
propKin[p_] := FourierDST[Ke*FourierDST[p, 1] , 1] ; 

(* propagation *) 

psil = propKin [propPot [Dt/ (2M) ,psi]] ; 

psi2 = Nest [propKin [propPot [Dt/M, #]] &, psil, M-l] ; 

propPot [Dt/ (2M) , psi2] ] 
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exercises 


Q4.7 Extend the split-step method of section 4.2 to generate not just the final state 
but all intermediate states as well. Hint: use the NestList command as in 
subsection 4.1.3 (page 85). 


Q4.8 Given a particle moving in the range xe [0,1] with the scaled non-linear Hamil- 
tonian 


- 1 d 2 

7r 2 d*2 +n 


r- 4\ 


-1 


+g|ig(x)| 2 , 


(4.43) 


W( x) 


do the following calculations: 


1. Plot the potential for Q = 1 and S - 4 (use g - 0). What are the main 
characteristics of this potential? Hint: compute V (4), 17(4 + 5), 1 7'(4±5). 

2. Calculate and plot the time-dependent density |t//(x, t)\ 2 for Q = 50 and 

g = 0, starting from i//o(x) ex exp [- with xq - 0.2694 and cr = 

0.0554. Calculate the probabilities for finding the particle in the left half 
(x < 4) and in the right half (x > 4) up to t — 100. What do you observe? 

3. What do you observe for n = 50 and g = 0.1? Why? 


4.2.1 imaginary- time propagation for finding the ground state of 
the non-linear Schrodinger equation 

You may remember from statistical mechanics that at temperature T, the density 
matrix of a system governed by a Hamiltonian is 


p{p) = Z~ l ip)e~^ (4.44) 

with p = l/(fc B T) in terms of the Boltzmann constant = 1.3806488(13) x 10 -23 J/K. 

The partition function Z(/l) = makes sure that the density matrix has the 

correct norm, Trp = 1. 

We know that at zero temperature the system will be in its ground state |y), 3 


lim pQ8) = |y>(y|. (4.45) 

/3— oo 


If we multiply this equation with an arbitrary state |i//> from the right, and assume 
that (y \y/) ^ 0, we find 

lim p(/3)|i//> = |y>(y|ig>. (4.46) 

(3—>oo 

The ground state is therefore 


ly> = 


limp^ 00 p(/3)[t^) 

<y|i//> 


1 

<yl y)Z{p) 


x lim e P^\y/). 

oo 


(4.47) 


This means that if we take (almost) any state \y/) and calculate limp^oo e we 

find a state that is proportional to the ground state. But we already know how to do 

3 For simplicity we assume here that the ground state is non-degenerate. 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 


88 


CHAPTER 4. REAL-SPACE SYSTEMS 


this: the wavefunction e~^\y/) is calculated from i//> by imaginary-time propaga- 
tion. In fact the algorithm of subsection 4.1.3 remains valid if we replace i{t-to)IH<-~ 
(’>, and so does its extension to the non-linear Schrodinger equation (section 4.2). 
The only caveat is that, while regular time propagation (subsection 4.1.3) is uni- 
tary, imaginary-time propagation is not. The wavefunction must therefore be re- 
normalized after each imaginary-time evolution step (with the Normalize func- 
tion), particularly before calculating the non-linear potential in the Gross-Pitaevskii 
equation. 

For a computer implementation we modify Equation (4.47) to 


jim e -M6p*\ )= am \ e -ep* ] M ly/) 

M—*oo M—oo L J 


Trotter Equation (4.31) 

i 


U P jp 

n £? 2 ^pot 


M 


lim \e 2 -^pote s P-^r<n e f 

M—*oo l J 

q F^POt '^ c in ^pot ... 

s. y > 


W) 


< 5/3 y7/> 

p 2 e ^ c P ( 


(M - 1) repetitions < 


for a fixed “imaginary-time” step 8fS, and iterate until the term in the square bracket 
no longer changes and the infinite-/! limit (M — ► oo) has effectively been reached. 


in [306] :=groundstate [db_?NumericQ , g_?NumericQ , 
tolerance_ : 10~ (-10) ] := 

Module [{Ke, psiO, psil ,psi2 ,propKin,propPot , gamma}, 

(* compute exponentials *) 

Ke = Exp[-db*Table[n~2,{n, l,nmax}]] ; 

(* propagate with the potential operator *) 
propPot [ddb_ ,p_] := 

Normalize [Exp [-ddb* (Wval+g* (nmax+l)*Abs [p] ~2) ] * p] ; 
(* propagate with the kinetic operator *) 
propKin[p_] := 

Normalize [FourierDST [Ke*Four ierDST [p , 1] , 1] ] ; 

(* random starting point *) 

psiO = Normalize @ RandomComplex[{-l-1 , 1+1} ,nmax] ; 

(* propagation *) 

psil = propKin [propPot [db/2, Normalize [psiO] ]] ; 
psi2 = FixedPoint [propKin [propPot [db,#]]&, psil, 

SameTest->Function [{pi ,p2} ,Norm[pl-p2] ctolerance] ] ; 
(* ground state *) 
gamma = propPot [db/2, psi2] ; 
gamma] 


The last argument, tolerance, is optional and is given the value 10 -10 if not spec- 
ified. The FixedPoint function is used to apply the imaginary-time propagation 
until the result no longer changes (two consecutive results are considered equal if 
the function given as SameTest returns a true result when applied to these two re- 
sults) . 
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The ground state of the time-independent non-linear Schrodinger equation sat- 
isfies 


2m dx 2 


+ V(x) + g \y/(x)\ 2 


V'(x) = py/ix), 


(4.49) 


where /i is called the chemical potential and takes the place of the ground-state 
energy in the time-independent linear Schrodinger equation. Integrating Equa- 
tion (4.49) by y/* (x) from the left and integrating over x gives 


'=/ 


i fr*{x) 


h 2 d 2 

- 2~ dx 2 + V M + ' 


It 

2 m 


/ 


dig(x) 


dx 


i//(x)dx 

dx + J [ V(x) + g|i^(x)| 2 ] |i//(x)| 2 dx, (4.50) 


where we have assumed that / [i//(x)[ 2 dx = 1. We use this to calculate the chemical 
potential in In [306] by replacing line 20 with 

(* chemical potential *) 

mu = Table [n~2, {n, 1 ,nmax}] . Abs [FourierDST [gamma, 1 ] ] ~2 
+ (Wval+g* (nmax+1) *Abs [gamma] ~2) . Abs [gamma] ~2 ; 

(* return ground state and chemical potential *) 

{mu, gamma}] 


and adding the local variable mu in line 3. 

exercises 

Q4.9 Given a particle moving in the range xe [0,1] with the scaled non-linear Hamil- 
tonian 

Id 2 [l ] 2 2 

-j-^r + 500 1 x - - I +g |1//(X)| 2 , (4.51) 

do the following calculations: 

1. For g - 0 calculate the exact ground state |f) (assuming that the parti- 
cle can move in x e 1) and its energy eigenvalue. Hint: assume ((x) = 

/ and find the value of a which minimizes 

2. Calculate the ground state lim^oo e~^\() and its chemical potential 
by imaginary-time propagation (with normalization of the wavefunction 
after each propagation step), using the code given above. 

3. Plot the ground-state wavefunction for different values of g. 

4. Plot the chemical potential as a function of g. 

4.3 several particles in one dimension: interactions 

We have seen in subsection 2.4.2 (page 40) how to describe quantum-mechanical 
systems with more than one degree of freedom. This method can be used for de- 
scribing several particles moving in one dimension. In the following we look at two 
examples of interacting particles. 



90 


CHAPTER 4. REAL-SPACE SYSTEMS 


4.3. 1 two particles in one dimension with contact interaction 


We first look at two particles moving in a one-dimensional square well of width a 
and interacting through a contact potential 8(x\ - x 2 ). Such potentials are a good 
approximation of the interactions taking place in cold dilute gases. The Hamiltonian 
of this system is 


Je= 


h 2 


2 m 


d 2 d 2 
dx 2 dx 2 

v ' 

'^kin 


+ HUi)+ H(x 2 ) + gd(xi-x 2 ), 


^pot 


^int 


(4.52) 


where l/(x) is the single-particle potential (as in section 4.1) and g is the interaction 
strength, often related to the s-wave scattering length a s . 

We describe this system with the tensor-product basis constructed from two 
finite-resolution position basis sets: 


\j 1 J 2 ) = \ji) ® 1 72> f° r 7 i»72 g {1,2,3, ...,rc max }. (4.53) 


Most of the matrix representations of the terms in Equation (4.52) are constructed 
as tensor products of the matrix representations of the corresponding single-particle 
representations since ^kin = ^kin,i®l + l®^kin ,2 andJ&p 0t = ^ po t,i«il + l®^pot, 2 - 
The only new element is the interaction Hamiltonian -7C\ m . Remembering that its 
formal operator definition is 

= [|xi>®l^2>]5(xi-x 2 )[<xi|«><x2l]dx: 1 dr2 (4.54) 

(while Equation (4.52) is merely a shorthand notation), we calculate its matrix ele- 
ments as 


0'l.72l^mtljl.J2> = S / 0'llM>(j2|X 2 )5(Xi -X2)<Xil7(>(X 2 |j2>dXidX2 

Jo 

~sj dj\ ix)dj 2 {x)dj^ {x)d jAxjdx. (4.55) 

We could in principle evaluate these integrals exactly, but notice that there are ^(^max) 
integrals to be computed, which quickly becomes unmanageably slow as « max grows. 
Instead, we can again do an approximation: since the basis functions Qj (x) are zero 
on all grid points except at Xj [see Equation (4.10)], the integrand in Equation (4.55) 
vanishes on all grid points xi, x 2 , . . ., x„ max unless j\ = j 2 = j\ — j' 2 - We thus approxi- 
mate the integral by zero if the j -values are not all equal: 


( J 1 1 ]2 1 >^mt I ]\ * j 2 ) 


■8: 


h.hJi.i'2 X % 


f 


■8j t (x)dx. 


(4.56) 


These integrals are not easy to do in all generality. Exact integration of the case j = 
nm If +1 , which is localized at the center of the square well (x « \ , if n max is odd) gives 


f Vim ax+l«dx=^ 
Jo 2 oa 


2(^max + 1) + ' 


1 + 1 


(4.57) 


We will use this expression to approximate all quartic overlap integrals Jq ( x)dx . 

2 fc2 

Mathematica code: we assume a = 1, express energies in units of E\ - * , and 

assume that the potential function W (x) is defined. First we define the grid size and 
the unit operator id acting on a single particle: 
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in [ 307 ] :=nmax = 10; 

in [308] :=xval = Range [nmax] / (nmax+1 ) ; 

in[309] : = id = SparseArray [Band [{1 , 1}] -> 1, { nmax, nmax }-] ; 


The total kinetic Hamiltonian is assembled via a Kronecker product (tensor product) 
of the two single-particle kinetic Hamiltonians: 

in[3io] :=HkinlM = SparseArray [Band [{1 , 1}] ->Range [nmax] ~2] ; 
in[3ii] :=HkinlP = FourierDST[HkinlM, 1]; 
in[3i2] :=HkinP = KroneckerProduct [HkinlP, id] 

+ KroneckerProduct [id, HkinlP]; 


The same for the potential Hamiltonian (here we assume no potential, that is, a 
square well): 

In[313] : = W [x_] = 0; 

In[314] : = Wval = W /@ xval ; 

in[3i5] :=HpotlP = Spar seArray [Band [{ 1 , 1 }] ->Wval] ; 
in[3i6] :=HpotP = KroneckerProduct [Hpot IP , id] 

+ KroneckerProduct [id, Hpot IP] ; 


The interaction Hamiltonian is only nonzero when j\ = j 2 = j[ = j' 2 , which can be 
represented with a SparseArray [{j_, j_, j_, j_}->l, {nmax, nmax, nmax, 
nmax}] massaged into the correct form: 

in[3i7] :=HintP = (2 (nmax+1 ) + 1/ (nmax+1 )) /3 * 

SparseArray [Flatten /<§ 

Flatten [SparseArray [{j_, j_, j_, j_} -> 1, 

{nmax, nmax, nmax, nmax}], 1]] 


The full Hamiltonian, in which the amplitude of the potential can be adjusted with 
the prefactor O, is 

in[3i8] :=HP [0mega_ , g_] = HkinP + 0mega*HpotP + g*HintP; 


We calculate eigenstates (the ground state, for example) with the methods already 
described previously. The resulting wavefunctions are in the tensor-product basis of 
Equation (4.53), and their corresponding densities can be plotted with 

in[3i9] :=plotdensity [r_] := Module [{rl ,r2}, 

(* make square array of density values *) 
rl = Part it ion [r, nmax]; 

(* add zeros at the limits *) 

r2 = SparseArray[{},{nmax+2,nmax+2}] ; 

r2 [ [2 ;; nmax+1 , 2;;nmax+l]] = rl; 

(* plot *) 

ListDensityPlot [r2 , DataRange -> {{0, 1}, {0, 1}}]] 


We thus plot a given wavefunction psi with 
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i in [320] := plot density [ (nmax+ 1 ) * Abs [psi] ~ 2 ] 


Here we plot the ground-state density for 0 = 0 (no potential, the particles move in 
a simple infinite square well) and g - +5 (repulsive interaction) as well as g = -2 
(attractive interaction), using n max =10 grid points: 


0=0, g=+5 



0.0 0.2 0.4 0.6 0.8 1.0 


0=0, g=-2 



*1 


*1 


We can see that for g > 0 the particles avoid each other, i.e., the density p{x\ , x 2 ) 
vanishes whenever x\ = X 2 , whereas for g < 0 they attract each other, i.e., the density 
p(x 1 , X 2 ) is maximal whenever x\ = X2. 

exercises 

Q4.10 Calculate the expectation value of the inter-particle distance, {x\ - X2), and its 
variance, {{x\ - X2) 2 ) - (x\ - X2) 2 , in the ground state as a function of g (still 
keeping Q = 0). Hint: The position operators xl and x2 are 


1 in [ 321 ] :=x = SparseArray [Band [{1 , 1}] ->xval] ; 

2 in [ 322 ] :=xl = KroneckerProduct [x , id]; 

3 in [ 323 ] :=x2 = KroneckerProduct [id , x] ; 


4.3.2 two particles in one dimension with arbitrary interaction 

Two particles in one dimension interacting via an arbitrary potential have a Hamil- 
tonian very similar to Equation (4.52), except that the interaction is now 

= V int {x\,x 2 ). (4.58) 

As an example, for the Coulomb interaction we have V; nt (xi, x 2 ) = r , | with Qi 

and Q 2 the electric charges of the two particles. For many realistic potentials 
only depends on |xi - X 2 I. 

The matrix elements of this interaction Hamiltonian can be approximated with 
a method similar to what we have already seen. The exact expression 

na 

<7l»j2l^intlj'i»72> = J o 5; 1 (Xl)5j 2 (X2)^in t (Xi,X2)5 J -;(Xi)^/(X 2 )dXidX2 (4.59) 
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is approximated by noticing that on the grid points, the numerator of the integrand 
is only nonzero if j\ = j\ and j 2 = j'/- 

0'l.72l^intl7nJ2> ~ 8 h,j[ 5 j 2 ,]' 2 X J 0 (Xl)d 2 j 2 (.X 2 Wmt(Xi,X 2 )dXidX 2 . (4.60) 

With dj(x) ~ 8{x-Xj) (in fact it is the best approximation to a Dirac d-function pos- 
sible in our finite-resolution basis), these matrix elements simplify to 

na 

Ui’ i' 2 ) ~ 8 ji,j[ 8 h,i' 2 x J o d(Xi - XjjSiXz - Xj 2 )V int {x 1 ,X 2 )dXidX 2 

= 8 h,j[ 8 h,j' 2 x f 7 ] nt ( -'C/ ] > ^/2 ) • (4.61) 

This is again very easy to evaluate without the need for integration over basis func- 
tions. But realistic interaction potentials are usually singular for x\ = x 2 (consider, 
for example, the Coulomb potential), and therefore this expression, Equation (4.61), 
fails for the evaluation of the matrix elements (j,j\^\at\j,j)- This problem cannot 
be solved in all generality, and we can either resort to more accurate integration 
(as in subsection 4.3.1) or we can replace the true interaction potential with a less 
singular version: for the Coulomb potential, we could for example use a truncated 
singularity for \x\ < A: 


Vmt{x) = 


Q1Q2 

4neQ 


x 


_j_ 

\x\ 

3A 2 — x- 
2A 3 


if |x| > A 
if \x\ < A 


(4.62) 


As long as the particles move at energies much smaller than Tj n i (+A) = ^5^ they 
cannot distinguish the true Coulomb potential from this truncated form. 


exercises 

Q4. 1 1 Consider two particles in an infinite square well, interacting via the truncated 
Coulomb potential of Equation (4.62). Calculate the expectation value of the 
inter-particle distance, {x\ - x 2 ), and its variance, {(x\ - x 2 ) 2 ) - <X| - x 2 ) 2 , in 
the ground state as a function of the Coulomb interaction strength (attractive 
and repulsive). Hint: set A = a/[n max + 1) in Equation (4.62). 


4.4 one particle in several dimensions 


An important application of the imaginary- time propagation method of subsection 4.2.1 
is the calculation of the shape of a three-dimensional Bose-Einstein condensate in 
a harmonic trap. In this section we use such a calculation as an example of how to 
extend single-particle lattice quantum mechanics to more spatial dimensions. 

The non-linear Hamiltonian describing a three-dimensional Bose-Einstein con- 
densate in a harmonic trap is 


J€=- 



( d 2 

2m 

\dx 2 


d 2 

+ dy 2 + 


d 2 ) 
dz 2 ) 


m 

H 

2 


oj 2 x x 2 + oj 2 y 2 + o) 2 z 2 j + (AT- 1) 


4 nh 2 i 


-| y/(x,y,z)\ 2 , 
(4.63) 


where we have assumed that the single-particle wavefunction y/{x,y,z) is normal- 
ized: f\y/{x,y,z)\ 2 dxdydz-l. 
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We perform this calculation in a square box, where \x\ < |, \y\ < |, and \z\ < 
we will need to choose a large enough so that the BEC fits into this box, but small 
enough so that we do not need an unreasonably large n max for the description of 
its wavefunction. Notice that this box is shifted by | compared to the [0 ... a] boxes 
used so far; this does not influence the calculations in any way. The energy scale of 
this box is E\ - ■ Introducing the dimensionless coordinates x = xl a, y = yl a, 

and z-z/a, and defining the dimensionless wavefunction ys{x,y,z) - a 3l2 ij/(x,y, z), 
the dimensionless ffamiltonian is found from Equation (4.63): 


~E~i 



/ d 2 d 2 d 2 1 
1 dx 2 dy 2 dz 2 J 


+ Cl 2 x 2 + El 2 y 2 + Cl 2 z 2 + [N - l)y\y/{x, y, z) | 2 , 


(4.64) 


where Q. x - etc. are the dimensionless trap frequencies and 7 = ^ is the 

dimensionless scattering length (interaction strength). 

The ground state of this dimensionless non-linear Hamiltonian of Equation (4.64) 
can be found by three-dimensional imaginary-time propagation, starting from (al- 
most) any arbitrary state. Here we assemble a Mathematica function groundstate 
which, given an initial state psiO and an imaginary time step db, propagates until 
the state is converged to the ground state. 

First we define the dimensionless parameters of the problem. We will be consid- 
ering N - 100 0 87 Rb atoms in a magnetic trap with trap frequencies co x - 2n x 1 15 Hz 
and a)y-oj z -2nx 540 Hz. The 87 Rb atoms are assumed to be in the \F = 1 ,Mp - 1) 
hyperhne ground state, where their .v-wave scattering length is a s - lOOAag (with 
ag = 52.9177pm the Bohr radius). 

in[324] :=With [{m = Quant ity [" 86 . 909187 u"] , 
a = Quantity["10 um"] , 
wx = 2 Pi Quantity ["115 Hz"], 
wy = 2 Pi Quantity ["540 Hz"], 
wz = 2 Pi Quantity ["540 Hz"], 
as = Quantity [" 100 . 4 a0"] , 

H = Quantity ["fi"] } , 

Ox = m wx a~2/(Pi H) ; 

Oy = m wy a~2/(Pi H) ; 

Oz = m wz a~2/(Pi H) ; 
g = 8 as/ (Pi a) ;] 


Next we define the grid on which the calculations will be done. In each Cartesian 
direction there are « max grid points icj — xval [ [ j ] ] : 

in [325] :=nmax = 41; 

in [326] :=xval = Range [nmax] / (nmax+1) - 1/2; 


We define the dimensionless harmonic-trap potential: the potential has its mini- 
mum at the center of the calculation box, i.e., atx=:y = .z = | . 

in[327]:=W[x_,y_,z_] = 0x~2*x~2 + 0y~2*y~2 + 0z~2*z~2; 


We only need the values of this potential on the grid points. To evaluate this, we 
build a three-dimensional array whose element Wval [ [ j x , j y , j z] ] is given by the 
grid-point value W [xval [ [ j x] ] , xval [ [ j y] ] , xval [ [ j z] ] ] : 
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is in[328] :=Wval = Table [W [xval [ [j x] ], xval [ [j y] ], xval [ [j z] ]] , 
i6 {jx,nmax}, {jy ,nmax]- , {jz.nmax}] ; 


We could also define this more efficiently through functional programming: 


17 


in[329] :=Wval = Outer [W, xval, xval, xval]; 


The structure of the three-dimensional Wval array of potential values mirrors the 
structure of the wavefunction that we will be using: any wavefunction psi will be a 
n max x n m ax x Umax array of coefficients in our finite-resolution position basis: 

^max 

\fr(x,y,z)= Yx psi [ [jx, jy , jz] ]dj x (x)dj y (y)dj z (z). (4.65) 

jx,jy,jz= i 

From Equation (4.10) we find that on the three-dimensional grid points the wave- 
function takes the values 

ty(Xjx>Xj y <Xh) = («max+ D 3/2 psi[[jx, jy, jz]]. (4.66) 


The norm of a wavefunction is 

r Cl r 1 ^max 

/ |i//(x,y,z)| 2 dxdydz = / \y/(x, y,z)\ 2 dxdydz = £ |psi [ [jx , jy , jz] ] | 2 

Jo Jo h,jy,h = i 

= Norm [Flatten [psi] ] "2, (4.67) 


from which we define a wavefunction normalization function 


in[330] :=nn [psi_] := psi/Norm [Flatten [psi]] 


The ground state calculation then goes by imaginary-time propagation, with step 
size db corresponding to an evolution e~ db,:/2//:i per step. The calculation is done 
for N — n particles. Notice that the Four ierDST function can do multi-dimensional 
discrete sine transforms, and therefore the kinetic -energy propagator can still be 
evaluated very efficiently. The last argument, tolerance, is optional and is given 
the value 10 -6 if not specified. 
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in[33i] :=groundstate [n_?NumericQ , db_?NumericQ , 
tolerance_ : 10~ (-6) ] := 

Module [{Ke ,propKin,propPot ,psi0 ,psil ,psi2 , Gamma, muKin,muPot ,mulnt ,mu} , 
(* kinetic propagator in momentum basis *) 

Ke = Table [Exp [-db* (nx~2+ny~2+nz~2) ] , 

{nx,nmax}, {ny,nmax}-, {nz,nmax}] //N; 

(* kinetic propagator in position basis *) 
propKin [psi_] := FourierDST [Ke*FourierDST [psi , 1] , 1] ; 

(* random starting point *) 

psiO = nnQRandomComplex [{-1-1 , 1+1}, {nmax, nmax, nmax}] ; 

(* potential propagator in position basis *) 
propPot [b_?NumericQ , psi_] := 

Exp [-b* (Wval+g* (n-1) * (nmax+1) ~3*Abs [psi] ~2 ) ] *psi ; 
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(* first evolution step *) 

psil = nn [propKin [propPot [db/2,nn[psi0]]]] ; 

(* iterate evolution until wavefuction converges *) 
psi2 = FixedPoint [nn [propKin [propPot [db, #]]]&, psil, 
SameTest -> (Norm[Flatten[#l-#2] ] < tolerance &)]; 
(* one last half-iteration *) 

Gamma = nn [propPot [db/2 , psi2]]; 

(* chemical potential *) 

muKin = Flatten [Table [nx~2+ny~2+nz~2 , 

{nx , nmax} , {ny , nmax} , {nz , nmax}] ] . 

Flatten [Abs [FourierDST [Gamma, 1] ] "2] ; 
muPot = Flatten [Wval] .Flatten [Abs [Gamma] ~2] ; 
mulnt = g* (n-1) * (nmax+1) "3 * 

Total [Flatten [Abs [Gamma] ~4] ] ; 
mu = muKin+muPot+muInt ; 

(* return ground state and chemical potential *) 

{mu, Gamma}] 


As an example, we calculate the ground state for N — 1000 atoms and a time step of 
db = 0.001: 


in [ 332 ] : = fpg , p} = groundstate [1000 , 0.001]; 


The chemical potential is 


In [ 333 ] : = pg 
Out [ 333 ] = 228. 421 


From this result we can, for example, calculate the expectation values X = (x),Y- 
( y >, Z = < z >, XX = (x 2 ), YY = < y 2 ), ZZ = (z 2 ). We could define coordinate arrays as 

in [334] :=xc = Table [xval [ [j x] ] , {jx,nmax}, {jy,nmax}, {jz,nmax}]; 

in [335] :=yc = Table [xval [ [j y] ] , {jx,nmax}, {jy,nmax}, {jz,nmax}]; 

in [ 336 ] : = zc = Table [xval [ [j z] ] , {jx,nmax}, {jy,nmax}, {jz,nmax}]; 


or we could define them more efficiently as follows: 


in [ 337 ] : = ones = Array [1&, 

nmax] ; 



in [ 338 ] :=xc = Outer [Times 

xval , 

ones , 

ones] ; 

in [339] :=yc = Outer [Times 

ones , 

xval , 

ones] ; 

in [340] : = zc = Outer [Times 

ones , 

ones , 

xval] ; 


The desired expectation values are then computed with 

in[ 34 i] :=X = Total [Flatten [xc * Abs[p]~2]]; 
in [342] :=Y = Total [Flatten [yc * Abs[p]~2]]; 
in [343] : = Z = Total [Flatten [zc * Abs[p]~2]]; 
in [344] :=XX = Total [Flatten [xc~2 * Abs[p]~2]]; 
in [345] :=YY = Total [Flatten [yc~2 * Abs[p]~2]]; 
in [ 346 ] : = ZZ = Total [Flatten [zc~2 * Abs[p]~2]]; 
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The size of the BEC is then calculated from these as the standard deviations of the 
position in the three Cartesian directions: 

, in[347] :={Sqrt [XX-X~2] , Sqrt [YY-Y~2] , Sqrt [ZZ-Z~2] > 

2 out [347] = {0. 158723, 0.0419921, 0. 0419921} 


exercises 


Q4.12 Take the BEC Hamiltonian of Equation (4.63) in the absence of interactions 
( a s - 0) and calculate analytically the expectation values (x 2 ), (y 2 ), (z 2 ) in the 
ground state. 


Q4.13 Take the BEC Hamiltonian of Equation (4.63) in the limit of strong interactions 
(Thomas-Fermi limit) where the kinetic energy can be neglected. The Gross- 
Pitaevskii equation is then 


y \o> 2 x x 2 + (o 2 y 2 + o)^z 2 j + (N - 1) 
which has two solutions: 


4nh 2 a s 

m 


\y/(x,y,z) | 2 


\fr[x, y, z) = py/{x, y, z), 
(4.68) 


Oor 

\y{x,y,z ) = j n-f(u, 2 x x?+(o 2 y y 2 +(olz 2 ) (4.69) 

Together with the conditions that \y/(x, y, z)\ 2 > 0, that i//(x, y, z) should be con- 
tinuous, and that f \y/(x, y, z)\ 2 dxdydz = 1, this gives us the Thomas-Fermi 
“inverted parabola” density 


P o 
0 


I y/[x,y,z)\‘ 
with the central density 


(tf-(i)-(tri *(*r + (ir *(*)*>■ 

if not, 


Po 


87 1 


225 m 6 u> 2 ci) 2 a) 2 
h 6 a 2 (N- l) 3 


(4.70) 


(4.71) 


the Thomas-Fermi radii 


15h 2 a s (N- 1 )ci) y (o z 


Rx — 


and the chemical potential 


R y = 


\Sh 2 a s {N - 1) a) z (i) x 

m 2 o)y 


, Rz = 


15 h 2 a s {N- 1) (m x (i)x 


(4.72) 


p - - 1 225 mh 4 a 2 {N - l) 2 a> 2 co 2 colj 5 . 


(4.73) 


Usingthis density, calculate the expectation values (x 2 ), (y 2 ), (z 2 ) in the ground 
state of the Thomas-Fermi approximation. 


98 


CHAPTER 4. REAL-SPACE SYSTEMS 


Q4.14 Compare the numerical expectation values < x 2 ), (y 2 ), (z 2 ) of our Mathematica 
code to the analytic results of Q4.12 and Q4.13. What is the maximum 87 Rb 
atom number N which allows a reasonably good description (in this specific 
trap) with the non-interacting solution? What is the minimum atom number 
which allows a reasonably good description with the Thomas-Fermi solution? 

Q4.15 Consider a 87 Rb Bose-Einstein condensate in a harmonic trap, described by 
the non-linear Hamiltonian of Equation (4.63). Take ojy-oj z -2nx 500 Hz 
and a scattering length a s = 100.4ao- Compute the expectation values ( x 2 ), 
( y 2 ), (z 2 ) for several values of oi x and try to interpret the asymptotes o> x 0 
and <i) x — > oo. 


Chapter 5 


combining space and spin 


In this chapter we put many of the techniques studied so far together: spin degrees 
of freedom (chapter 3) and spatial degrees of freedom (chapter 4) are combined with 
the tensor-product formalism (chapter 2). 


5.1 one particle with spin in one dimension 

5.1.1 separable Hamiltonian 

The simplest problem combining a spatial and a spin degree of freedom in a mean- 
ingful way consists of a single spin- 1/2 particle moving in one dimension in a state- 
selective potential: 


Je= 


h 2 d 2 
2 m dx 2 


+ Vo(x) + V z W §z. 


(5.1) 


where S z = |<r z is given by the Pauli matrix. As was said before, Equation (5.1) is a 
short-hand notation of the full Hamiltonian 

u2 noO noO rOO 

— / dx|x) — 5 -(x|®l+ / dx|x) Vo(x)(x| ® 1 + / dx|x) V z (x)(x\ ® S z , 

2 m J—oo dx z J - oo J- oo 

(5.2) 

where it is more evident that the first two terms act only on the spatial part of the 
wavefunction, while the third term couples the two degrees of freedom. 

The Hilbert space of this particle consists of a one-dimensional degree of free- 
dom x, which we had described in chapter 4 with a basis built from square-well 
eigenstates, and a spin- 1/2 degree of freedom S — described in the Dicke basis 
(chapter 3). This tensor-product structure of the Hilbert space allows us to simplify 
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the matrix elements of the Hamiltonian by factoring out the spin degree of freedom, 

h 2 


u r°° r°° 1 r°° 

<<M T> = -— / rU)/'Wdx<m>+ 0*U)Vb(*)VWdx<ili> + - / 0*(x)V z (x)i/dx)dx(tl<7 z lt> 

ZJTl J — oo J — oo ^ J—oo 

n 2 


X 


—oo 

■oo 


2m J- o 
H 2 '■°° 


■/. 


c p (x)i/a (x)dx + J <p [x)Vo(x)y/[x)dx + 


1 r c 

2 J-c 


<p (x)V z {x)xf/(x)dx 


2 m 


f 


■L 


i 


(<M l^flVh l> = - — I (f> ix)y/ (x)dx(f ||) + / (p (x)V r 0 (x)i//(x)dx(tli> + - / <p (x)V z (x)t^(xldx<||d- z li> 


f 

2 J-c 


= 0 


<0 


fe2 r oo ooo i roo 

,\\je\y,\) = -— 0*(x)i//'(x)dx<||t> + / ^*(x)Vo(x)i//(x)dx<|||> + - / 0* (x) V z (x)i^(x)dx<||(j z ||> 

Z/?Z J—oo J—oo Z J—oo 


= 0 


fi 2 


r 

J— c 


2m 
fi 2 
2 m J-, 


X 


((f), l \J?\y/, !> = -— I (p My/ (x)dx<|||> + / cp (x)V 0 (x)y/(x)dx(l\l) + - / cp (x)V z (x)i//-(x)dx<||<7 z ||> 


If 

2j- 0 


/ oo r*oo | r*co 

0*(x)i// , (x)dx+ / (/>* (x) Vo(x)i^(x)dx — / </>* (x) V z (x)i//'(x)dx. 

-00 J—oo 2 J—oo 


(5.3) 


We see that this Hamiltonian does not mix states with different spin states (since 
all matrix elements where the spin state differs between the left and right side are 
equal to zero). We can therefore solve the two disconnected problems of finding the 
particle’s behavior with spin up or with spin down, with effective Hamiltonians 


H 2 d 2 1 

= _ 2m dx^ + Vq{x) + -V z (x), 

H 2 d 2 1 

^-^ + VoM -2 VM - 


(5.4) a 

(5.4) b 


These Hamiltonians now only describe the spatial degree of freedom, and the meth- 
ods of chapter 4 can be used without further modifications. 


5.1.2 non-separable Hamiltonian 

A more interesting situation arises when the Hamiltonian is not separable as in sub- 
section 5.1.1. Take, for example, the Hamiltonian of Equation (5.1) in the presence 
of a transverse magnetic field B x , 

h 2 d 2 - 

■?€= + v o( x] + v z(x)S z + B X S X . (5.5) 

The interaction Hamiltonian with the magnetic field is not separable: 


((p,UB x S x \V,}) = ^B x 
((p,UB x S x \V,\) = ^B x 
(<P> I \B x S x \y/, f> = ^B x 
((p,l\B x S x ty,\) = ^B x 



0*(x)l//(x)dx(f |<T;c[f> = 0 

1 f°° 

0*(x)i//(x)dx(f|<7;c[|> = -B x I <p* (x)y/(x)dx 

^ J—OO 

1 r°° 

0*(x)i//(x)dx(||<7;c[f> = -B x I <p* (x)y/(x)dx 

^ J—OO 

0*(x)i//(x)dx(i|(T x |i) = 0. 


(5.6) 
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Therefore we can no longer study separate Hamiltonians as in Equation (5.4), and 
we must instead study the joint system of spatial motion and spin. In what follows 
we study a simple example of such a Hamiltonian, both analytically and numerically. 
We take the trapping potential to be harmonic, 


VoM = -rmo 2 x 2 


(5.7) 


and the state-selective potential as a homogeneous force, 


V z {x) = ~Fx. 


(5.8) 


ground state for B x - 0 

For B x - 0 we know that the ground states of the two spin sectors are the ground 
states of the effective Hamiltonians of Equation (5.4), which are Gaussians: 

( x-fi ,2 

e l 2a I 

(x\y\) = — =^= ® It) <x[yj> = — (5.9) 

\J asfTn v <j\/2n 


with p = 2 maP an3 <T ~ V 2 mw • These two ground states are degenerate, with energy 

E— ^htj) - ^ 2 . In both of these ground states the spatial and spin degrees of free- 
dom are entangled: the particle is more likely to be detected in the ||) state on the 
right side (x > 0) , and more likely to be detected in the [ j ) state on the left side (x < 0) 
of the trap. This results in a positive expectation value of the operator x ® S z : 

<yj|x® S z |yj> = <n|x® S-|yj> = | = rc 2 . (5.10) 


perturbative ground state for B x > 0 

For small | B x | the ground state can be described by a linear combination of the states 
in Equation (5.9). If we set 


[yp> = ax|y T > + ^x|yj> (5.11) 

with | a | 2 + | p | 2 = 1 , we find that the expectation value of the energy is 


<7pl^lrp> = ial 2 <rt \^\r\) + a*P(r\ + P* a ( n l^l7t> + \P\ 2 (n l^ln) 

1 F 2 1 f 2 

= -hxo T + -B x (a*B+l3*a)e4mh^ (5.12) 

2 8 mio z 2 


For B x > 0 this energy is minimized for a - 1 / x/2 and /) = - lly/2, and the perturba- 
tive ground state is therefore the anti-symmetric combination of the states in Equa- 
tion (5.9) 

r*-n 2 -(ihl) 2 

(x|y p ) = ——= s IT) - -= <8>ll). (5.13) 

\2a\/2n \/2o\j2n 


with energy 


<7pl^lrp> = \ hw ~ 


F 2 

8 raw 2 


— By6 4 m / mo 3 

2 


(5.14) 
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The energy splitting between this ground state and the first excited state, 


-lizE ) 1 2 

e i 2 a I e 

<*|e p > = ; , <g> I T) + 


3 ? 


\/ 2<j\/2n \J 2o\/2n 


®ll>- 


(5.15) 


is A E — (c p |^f|e p ) - < 7 p[^if[ 7 p) = B x e unhw 3 , which can be very small for large expo- 


nents 


4 mhiD 3 ' 


numerical calculation of the ground state 


For a numerical description of this particle we first re-scale the Hamiltonian to elim- 
inate unnecessary units. As usual we describe the spatial degree of freedom in a box 
of size a, with energy scale E\ = * ; we set x- ax and use the range -| < x < |.' 

The scaled Hamiltonian is 




1 d 2 
jt 2 dx 2 


+ Q 2 x 2 - fx ® S z + b x S x 


(5.16) 


with Q = a )'- and b x - B x the dimensionless parameters of the 
problem. 

We describe the spatial degree of freedom with the finite-resolution position ba- 
sis of section 4.1.1: 


in [348] :=nmax = 20; 

in [ 349 ] :=xva.l = Range [nmax] / (nmax+1 )- 1 /2 ; 


The operator x is approximately diagonal in this representation: 


in[350] :=xop = Spar seArray [Band [{ 1 , 1 }] -> xval] ; 


The identity operator on the spatial degree of freedom is 


in [ 351 ] : = idx = SparseArray [Band [{1 , 1}] -> 1, {nmax, nmax} ] ; 


The Pauli operators for the spin degree of freedom are 


In [352] : = ids = {{1 , 0} , {0 , 1}} ; 

in [ 353 ] : = {sx, sy , sz}=Table [SparseArray [PauliMatrix [i] / 2] , {i , 1 , 3}] ; 


The kinetic energy operator is constructed via a discrete sine transform, as before: 


in [ 354 ] :=HkinM = SparseArray [Band [{1,1}] ->Range [nmax] ~2] ; 
in [ 355 ] :=HkinP = FourierDST [HkinM , 1] ; 


From these we assemble the Hamiltonian: 

1 Until now we had always used 0 < x < 1. Shifting this domain to — \ < x < | does not change anything 

in the computational methods presented so far. 
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in[356] :=H [0mega_ , f_, bx_] = 

in [357] := KroneckerProduct [HkinP , ids] 

in [358] := + 0mega~2 * KroneckerProduct [xop . xop , ids] 

in [359] := - f * KroneckerProduct [xop , sz] 

in [360] := + bx * KroneckerProduct [idx , sx] ; 


We compute the ground state of this Hamiltonian with 


in[36i] : = Clear [gs] ; 

in[362] :=gs [Omega_?NumericQ, f_?NumericQ, bx_?NumericQ] : = 
gs [Omega, f, bx] = 

-Eigensystem [-H [N [Omega] , N [f ] , N [bx] ] , 1 , 

Method -> {"Arnoldi", "Criteria" -> "RealPart" , 
Maxlterations -> 10~6}] 


Once a ground state \ j) has been calculated, for example with 

in [363] := gamma = gs[100, 10000, 1000] [[2, 1]]; 
in [364] :=Dimensions [gamma] 

Out [364] = {40} 


the usual problem arises of how to display and interpret the wavefunction. In or- 
der to facilitate this analysis, we first re-shape the ground state to better reflect the 
tensor-product structure of our Hilbert space: 

in [365] :=gammaA = Partition [gamma, 2]; 
in [366] :=Dimensions [gammaA] 

Out [366] = {20, 2} 


In this way, gammaA [ [j , s] ] is the coefficient corresponding to the basis function 
| j) ® || - s) = \ j, | - s), with the definitions | + |> = [|> and | - |) = |[) for the spin 
part (so that 5=1 corresponds to the |f) and s = 2 to the [[) state). From this re- 
shaped ground state we calculate the re-shaped density matrix 

in[367] :=rhoA = 0uter[Times, gammaA, Conjugate [gammaA] ] ; 
in [368] :=Dimensions [rhoA] 
out [368] = {20, 2, 20, 2} 


In this density matrix, Cj liSlt j 2iS2 = rhoA [ [j 1 , si , j 2 , s2] ] is the coefficient corre- 
sponding to the basis function \ji, | — Si)(j 2 , \ - .sy ( in the basis expansion of the 
density matrix: 


ftmax 2 Wmax 2 3 3 

P= E E E L C h ,s 1 ,j2,s 2 \h^-Sl)(j2,--S2\. 
ji=l S!=l j 2 =l s 2 =l z z 


(5.17) 


Specifically, (nmax+1) *rhoA [ [j , s , j , s] ] is the probability of detecting the par- 
ticle at Xj in spin state | - s. With this density matrix we calculate the following 
quantities: 
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Spin-specific densities: if we only detect particles in spin state ||>, we measure in 
effect the expectation values of the spin-selective density operator p\ (x,y) = 
\x){y\ ® |I)<I|. The spin-selective density values are found from the expansion 
of Equation (5.17) and the trace Tr(A) = ( j , § - s\A\j, | - s): 


p t (x, y) = Tr[p • pj (x, y)] 

n max 2 ftmax 2 3 3 

= Tr L L L L c h.si.h,s2\h’? - siX/2, - - s 2 | - |^)<y[ <s ITXTI 

Jl = lSl = l 72=1 s 2 =l Z Z 

^max 2 AZmax 2 AAmax 2 3 3 3 3 

= L L L L L L c ;i,sij 2 ,s 2 0'-^-sl7'i.^-si>0'2,--S2|x ) t><y ) T l7,--s> 

i=i s =i h= ls i =1 12=1 S 2 =i z z z z 

^max ^max 

- £ £ c h , lk yd j2 {x)d h {y). (5.18) 

71 = 1 ;2=1 


Specifically, if x — Xj x and y = xj y lie exactly on grid points of our finite-resolution 
calculation grid, then from Equation (4.10) that 

P\iXj x ,Xj y ) = («max+ l)Cj x ,l,j y ,\. (5.19) 

That is, the detected density of spin-up particles is (at least on the grid points) 
given directly by the coefficients of rhoA computed above: 


i in[369] :=rhoup = (nmax+1) * rhoA[[All, 1, All, 1]]; 


In the same way the density of particles in the spin state | J. ) is 


i in[370] :=rhodown = (nmax+1) * rhoA[[All, 2, All, 2]]; 


They are plotted with a similar function to what was used on page 91, 


1 in[37i] :=plotdensity [r_] := Module [{r2} , 

2 (* add zeros at the limits *) 

3 r2 = SparseArray [-Q , {nmax+2 , nmax+2}] ; 

4 r2 [ [2; ; nmax+1 , 2;;nmax+l]] = r; 

5 (* plot *) 

e ListDensityPlot [r2, 

7 DataRange -> {{-1/2, 1/2} , {-1/2, l/2»]] 


1 in[372] :=plotdensity [rhoup] 

2 in[373] :=plotdensity [rhodown] 
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Pt(x.x') Pj.(x,x') 



- 0.4 - 0.2 0.0 0.2 0.4 - 0.4 - 0.2 0.0 0.2 0.4 


x x 

Reduced density matrix of the spatial degree of freedom: using Equation (3.36) we 
“trace out” the spin degree of freedom to find the density matrix in the spatial 
coordinate: 

i in [374] : = rhox = (nmax+1) * Sum[rhoA[[All,s,All,s]] , {s,l,2}]; 


This density is just the sum of the spin-specific densities shown above. 

P(x.x’) 



- 0.4 - 0.2 0.0 0.2 0.4 


x 


Reduced density matrix of the spin degree of freedom: we can do the same for the 
reduced matrix of the spin degree of freedom: 


1 in[375] :=rhos = Sum [rhoA [ [ j , All, j, All]], {j , 1, nmax}] 

2 Out [375] = {{0.5, -0 . 205979} , {-0.205979, 0.5}} 


Spin expectation value: if we only detect particles at a given position x, the expec- 
tation value for the measured spin is given by (<t z (x)> = + : 


1 in [376] : = avgs = Table [Sum [rhoA [ [j , s , j , s] ] * (3/2-s) , {s,l,2}]/ 

2 Sum [rhoA [ [j ,s,j ,s]] , {s,l,2}], {j,l,nmax}]; 

3 in[377] :=ListLinePlot [avgs , PlotRange -> All, 

4 DataRange -> {xval[[l]], xval[[-l]]}] 
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This graph confirms the observation that particles detected on the left side are 
more likely to be in the [|> state, while particles detected on the right side are 
more likely to be in the [f) state. 

5.1.3 exercises 

Q5.1 In the problem described by the Hamiltonian of Equation (5.5), calculate the 
following expectation values (numerically) for several parameter sets {Q, /, b f \: 

• ( x ) for particles detected in the | f ) state 

• ( x ) for particles detected in the 1 1) state 

• ( x ) for particles detected in any spin state 

• the mean and variance of x ® S z 


Chapter 6 

path-integral methods 


With the approximations associated with the finite-resolution position basis set (sec- 
tion 4.1.1) we have significantly reduced the complexity of performing calculations 
of quantum-mechanical systems with continuous degrees of freedom such as their 
motion in space. When we study a very large number of particles, however, these ap- 
proximations are still not sufficient and computers are still overwhelmed. Consider, 
for example, the extension of the problem of subsection 4.3.1 to N particles moving 
in three-dimensional space. Representing any wavefunction of this system in the 
position-basis requires /r^ ax complex numbers; for N - 20 and /? max = 20, which are 
both not very large numbers, we already require about 10 78 complex numbers for 
the complete description, which approximates the number of particles in the uni- 
verse. 

The Trotter decomposition of Equation (4.31) we used in subsection 4.1.3 (real- 
time dynamics) and subsection 4.2.1 (imaginary-time dynamics) lends itself to a 
quantum-mechanical description that circumvents this problem and can be used 
to estimate expectation values without calculating the full wavefunction or density 
matrix. And since we are ultimately interested in expectation values (measurable 
quantities), not in wavefunctions and density matrices (unmeasurable representa- 
tions), this can be of significant use. 


6. 1 path integrals for propagation in time 

Assume for a moment that we study a system with a time -independent Hamilto- 
nian Equation (2.33) gives an explicit expression for the solution of the time- 
dependent Schrodinger equation through the propagator ^(f) = exp(-i J€tlK). 

Here we wish to calculate matrix elements of this propagator in the position ba- 
sis, 

{x!\e~^ tlh \x), (6.1) 

where both x, and x are configuration vectors describing the 3 A spatial coordinates 
of the system; x = {x\,y\,z\,X 2 ,y 2 < Z 2 , . . . , jcjv, yjv, Zjvb The matrix element of Equa- 
tion (6.1) computes the amplitude with which a state (configuration) x turns into a 
state (configuration) x' during an evolution time t. It will be useful to look at the 
points x and x as the starting and end points of a path through 3,'V-dimensional 
configuration space. 
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First we apply the Trotter expansion of Equation (4.31) to the propagator: 

V I M—oo 1 > 


1 M 


— if jp _ if dp _ it dp _ it dp _ if db _ if dp _ if dp 

— li m g 2Mh^P°t g Mh^king Mh^ pot . . . g Marking Mh^P 01 g Mh^^n g zMh^P' 
M— ™ s ' 


if dp _ _ if dp 

(M - 1) repetitions of e Mh kin e Mh P ( 


Now we insert the unit operator 


/' 


1 = / |x>(J| d 3N x 


r 


■/. 


|xi)<jci|cki ® lyiXyildyi® / |zi><ziidzi® 


f 


■f 


|Z]v><Z]vldzjv (6.3) 


between each two operators in Equation (6.2). This unit operator integrates over all 
coordinates of all particles (3 IV coordinates in total), and the vector x represents all 
of these 3N coordinates; [ x) - \x\) ® |yi) ® \z\) ® • • • ® [zjv). With these insertions, 
matrix elements of Equation (6.2) become 


(x'\e~ i ^’ t,h \x) = lim f A^p«|je M ) 

M—ooJ 


x(x M \e Mh-^v° l \x M _i) ■ ■ ■ (x 2 \e Mh-^^x'^ix'^e MS^P ot |xi> 

— if dp. . _ if 

(M- 1) repetitions of <^ m +i|e Mh ^\x m ){x m \e Mh P oi \x m ) for m - (M- 1)... 1 

X <jei|e“ra^l4)(4|e“W^pot|j 0 >d 3W Jod 3W jeid 3N jc;---d 3N je' M _ 1 d 3N je M , (6.4) 


where we have set 3to = x and x' M = x as the starting and end points of the path. 
Equation (6.4) contains two kinds of integration variables, x m and x' m , which can 
be set equal because the potential Hamiltonian (including interaction potentials) 
is usually diagonal in the position representation: ^ pot = / l/(x)|jc)(x|d :i,v .r, and 
therefore 

{x' m \e~^^ ol \x m ) =6(x m -x' m )e~™ V( - Xm \ (6.5) 

Inserting Equation (6.5) into Equation (6.4) gives 


(x'\e l ^ tlh \x)= lim / e 2 mh V( - X m ' 1 

x {x M \e~^^ ]dn \xM-i)e~™ v< ' XM - l) ■■■(x 2 \e~^'^’ kin \x 1 )e~Wi V{Xl) 

Si ^ 

(M- 1) repetitions of (x m+ i\e~ Mh'bin \x m )e~ Mb v ^m) for m = (Af- 1) ...1 

X (xi\e~^^^\x 0 )e~2m v ^ o) d 3N Xi---d 3N x M -i> (6-6) 

where xq - x is the starting point and xm - x is the end point of the path. A fur- 
ther simplification of Equation (6.6) comes from the exact evaluation of the kinetic- 
energy matrix elements. If we assume that the kinetic energy represents the free 
motion of N particles of masses m n in three dimensions, 

*L h 2 ( d 2 d 2 d 2 ) 

•^kin — — X o ~T~ 2 + d~2 + d~2 ’ 

n=l 2m n \dxi dyfi dz i> 


(6.7) 
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then the 3 N coordinates of the particles propagate independently under . There- 

fore we first evaluate the matrix element for a single degree of freedom. With the 
conversion between the position basis \x) and the momentum basis |/c) given by a 
Fourier transform, 

1 r°° . 1 r°° 

I x) = —— I dke~ lkx \k ) Ik)-— —I dxe lkx \x ) (6.8) 

y/2n J—oo \j2n J—oo 

we find for a single coordinate (with a = jjjjm ) 


(x'\e ia fa 2 \x) = (x'|e‘“d* 2 


!\ 


i 


\/2n 


I. 


dke~ ikx \k) = 


1 r 
\ f 2n J-r 


d/c e 


-iak 2 -ikx] ^ 


oo 


i r c 

V2n J - c 


dk 


i ct o — \kx 
e ex 2 e 


(x'\ k) 


V2 h 


L f 

',TT J- 


= —l dke~ iak2 e~ ikx e ikx ' = 


— 1 i[x-x r ) 


47ra 


e , 


(6.9) 


} a jhz 0 -ikx _ [ y-°° tin)" d 2n ] „-i kx _ t-oo (ia)” r : „-\kx _ y-oo (-iafc 2 ) n „-i kx 

e - [2-„ =0 „! ~ 2 -n=0 nl *■ 1K) e ~ ^n= 0 n! e 


since e 

£ -iak 2 e ~ikx Therefore the full kinetic -energy propagator is 


<x'\e-^«\x) = n ( 2 ex P ( im ' lM[{Xn ~ X ' n)2 + iyn ~ y ' n)2 + [Zn ~ ^ 


n = 1 


2 nht 


2ht 


( 6 . 10 ) 

In the case where all particles have the same mass ( m n - m\/n) this propagator can 
be simplified to 


3JV , „ 

if ™ _ ( -imM\ 2 [ imM\\x- x \\^ 

exp 1 2 nt 


( 6 . 11 ) 


Using this latter form, Equation (6.11), for simplicity (a generalization to particles 
of unequal masses is straightforward but more complex), the expectation value of 
Equation (6.6) becomes 


3NM 

(x'\e~ i ^‘ tlH \x)= lim 

M—*oo 

X e w 1,2 e - ™ v9m - i) • • • e W II s 2-Si II 2 e - 4 mi 

11 v ' 

(M-l) repetitions ofeW l '»‘+ 1_ "" |2 TBi l ' S '" 1 for m = (M- 1) ... 1 

X e w 1*1 » V ^*o) d 3iv 3£i . . . d 31 V x (6 12) 

Notice that in this form, the integrand does not contain any operators any more: it 
is simply a product of numbers. This means that we can gather them as a single 
exponential: 


-i mM\ 2 
2nht j 


J 


e 2 Mh V &m) 


(x'\e 


-iJftlh 


\x) = lim 


. 3 NM 

-i mM\ 2 


m—o o\ 2nht I 


I 


exp 


it 

Jdh 


l M-l i 

-V(Xq)+ + -V&m) 


imM M* _ 2 

+ 2^ \\ x m~ x m-l\\ 

zrit m=1 


m-l 


d 3 N Xi---d 3 N x M -i- (6.13) 
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Now we look at the sequence Xo , 3c i , . . . , Xm as a path through the space of configura- 
tions of our ;V- particle system. Since we are considering the limit M —>o o, this path 
will eventually become continuous. We let x(j) describe this continuous path, with 
r = mt / M and x [mtIM) — x m . The starting point is x(0) = x, and x(f) = x is the end 
point of our path. With this definition we can make the substitutions 


lim — 
M—oo M 


l M-l 1 

-V(x 0 ) + Y V &m) + -V&m) 

^ m= 1 ^ 



V[x(T)]dr 


(6.14) 


(trapezoidal integration) and 


m M r t • 

lim — Yj \\Xm-Xm- 1 II 2 = I P(T)|| 2 dr (6.15) 

M ~°° t m= i Jo 


(tacitly assuming that the path x(r) is differentiable). With these substitutions we 
re-write Equation (6.13) as 


(x'\e~ ijetin \x) 


= lim 

M—*oc 


3 NM 

-i mM\ 2 


2nh t 


| J exp | ^ £ [y ||x(t)|| 2 - V[*(t)]| 


dr 


d 6N x 1 •••d 3iv x M _ 1 . 

(6.16) 


We recognize the integrand in the exponential of Equation (6.16) as the Lagrangian, 


m 


T£{x,x)='-\m 2 -V{x), 
and its integral as the action of the given path x{-), 


(6.17) 


SP[x ( •)] = f J£(x(t),x(t))cIt. (6.18) 

Jo 

With this we find the final form of the matrix element of the propagator, 

. rx' ■ 

<J'[e -uraft [3c>= / eK 5 ‘ w,1 ® f [jc(-)]. (6.19) 

Jx 


The symbol f~ 3> t [x{-)] denotes an integral over all (continuous and differentiable) 
paths x(t) running through 3/V-dimensional configuration space from x(0) = x to 

3 ATM 

x[t) — x 1 . The pre-factor ( ) " of Equation (6.16) has been absorbed into this 
symbol, and in general path integrals as in Equation (6.19) can only be interpreted 
up to a constant proportionality factor. When we calculate expectation values of 
operators, this is of no concern, as the following applications will show. 

Consider now what we have achieved: starting from a quantum-mechanical ex- 
pression for a matrix element in Equation (6.1), we have used the Trotter expan- 
sion to arrive at a numerical integral, Equation (6.19), that makes no reference to 
quantum mechanics at all. The price we pay is that Equation (6.19) requires us to 
find all continuous and differentiable paths which take us from the initial config- 
uration x to the final configuration x within time t. In this sense, the quantum- 
mechanical propagation from one state to another goes through all possible paths 
simultaneously; the amplitudes of all these paths, given by their action integral, can 
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interfere, as described by Equation (6.19). This gives a very intuitive picture to ex- 
periments such as Young’s double-slit, where the question of which slit the photon 
passes through is answered straightforwardly: it passes through both slits, and the 
interference pattern results from the interference of the action integrals of the two 
paths. 

While the path integral formula of Equation (6.19) looks clean and simple, it 
is not at all clear how such an integration is to be done in practice. For practi- 
cal calculations, we return to Equation (6.13) and use a finite number M of “time 
slices”. But if at each time slice we must integrate over the system’s spatial 3N co- 
ordinates (the (T n x.i . . . d iN x M -i terms), this path integration is impossible for any 
reasonably-sized problem. So what have we gained? We have gained in that Equa- 
tion (6.13) is still a simple numerical integration (albeit with very many integration 
variables) and can therefore be approximated by powerful techniques for evaluat- 
ing high-dimensional definite integrals. A commonly used technique is a stochastic 
Monte-Carlo evaluation of the path integral (the “Path-Integral Monte-Carlo” tech- 
nique): instead of summing over all paths x(t) connecting the starting point with 
the end point, we try to randomly generate a representative sample of paths, for ex- 
ample with the Metropolis-Hastings algorithm. 


6.2 path integrals for propagation in imaginary time 

With the substitution t — -i/t/i = -\h! (k\>,T), as in subsection 4.2.1 (page 87), we can 
calculate matrix elements of the thermal density matrix from Equation (6.13), 


(x'\e P^\x) = lim 


v 3 NM 

1 mM } 2 

r 

P 

\2nh 2 p) 

exp 

“m( 


1 M - 1 x 

-^(*o) + Y V{x m ) + -V{x M ) 


mM 


M 




Y \\Xm~X 


m- 1 II 


m-l 


d 3N Xi---d il ' l XM~i- (6.20) 


j3 n 


We again interpret the sequence xq,xi,..., xm as a path through the space of config- 
urations of our ,'V-particle system. We let x(t) describe this continuous path, with 
r = mfi/M and x{mj3IM) - x m . The starting point is x(0) = x, and x(/f) = x is the 
end point of our path. With this definition we can make the substitutions 



M— OO M 


1 M-l X 

-V(X 0 )+ Y V &m) + -V{X M ) 

2 m = 1 1 


P 

V[x(T)]dT 


( 6 . 21 ) 


and 


giving 


m M rP . 

lim T Y ll^m-^m-lll 2 = I ll*(T)|| 2 dr, 
P Jo 


( 6 . 22 ) 


(x'\e~P^\x) 


= lim 


3 MM „ 

l mM \ 2 r \ fP I m i , ^ 

J exp rJo fe"*™ * v,xM< ) 


dr 


j 3 iV- h 3 N~ 

a xi-”d xjvr-i- 
(6.23) 


Notice the modified sign in the integrand, compared to Equation (6.16). 
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6.2. 1 example: a single particle in a 1 D harmonic oscillator 

As a first example we study the harmonic oscillator Hamiltonian 

H 2 d 2 1 22 

— » + -m.br xr. 

2m dx 2 2 


(6.24) 


In what follows we calculate thermal matrix elements (x'\e ^\x)\ the calculation 
of propagator matrix elements follows the same scheme. 

summation over exact eigenstates 

We can exactly diagonalize Equation (6.24) with J€\ n) = E n \n), where the energies 

a 

I I_J fyj -X 

are E n - ho)[n + ±) and the eigenfunctions are {x\ n) = e 2 * 2 ' n terms of the 

Hermite polynomials H n (z) and with the length scale x - \J . The exact thermal 
matrix elements are therefore 


(x'\e~ p ^\x)= £ (x'\n')(n'\e~^\n)(n\x) 

n,n '= 0 


1 


= Y {x'\n)e P E "(n\x ) = — — y 

to iy/nto 2 "' l! 

1 -^ e -\ 


H n {xlx)H n (x'lx) c -^ c -0Hw( n +b 


1-/3 


x\fn 


72=0 


n\ 


exp 


(a#) 2 tanh(|)-(jg;) 2 coth(|) 


Xy/2nsmh{() 


, (6.25) 


where we abbreviate ( — fihoj — jrj, as the scaled inverse temperature, and we have 
made use of the identity 


E 

72=0 


2 w (2 w {xr + yr ) -2xy 

H n {x)H n {y)w n e 


n\ 


V 1 — 4 w 2 


(6.26) 


The partition function is 

~ roo pc 

Z(/3) = Tr(e~^) = (x\e~ fije \x)dx = 

J — <30 J — C 


exp | - (f) 2 tanh(|)] 


X\/27rsinh(<n 


dx 


1 


= - csch ( - 

2 [2 


(6.27) 


in terms of csch(z) = 1 / sinh(z), and hence the density matrix is 

,\ z fr\ ( m 2 


<x'|p(/3)|x> = 


<x'|e P*\x) 

Tr(e-A^) 


exp 


■(¥) tanh(|)-(^) coth(l) 


7T COth ^ j 


(6.28) 


We will compare our path integral calculations with this exact result. 
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path integral formulation 

In order to express the thermal density matrix as a path integral, we return to Equa- 
tion (6.20): for our one-dimensional problem (3 N >->• 1), 


(x'\e~^\x)= lim 


I mM \ 2 


m — oo\2nh 2 p ) 


f 


exp 


raor/3 
2 M 


i M - 1 
1 2 . V- J2. 


^o 2 +l 


X" + -X 


m= 1 


M 


mM M. ? 


dxi---dx M -i. (6.29) 


Here is an example of three different concrete paths for M - 5, starting at x - xq - 
0.42 and ending at x! — x$ — 0.14, passing through four intermediate points (X] , x 2 , x 2 , x 4 ) 



0.0 0.2 0.4 0.6 0.8 1.0 


x(r) 

For M — 5 we would now have to integrate over all possible intermediate points 
[xi,X 2 ,X 3 ,Xi), performing a four-dimensional integral, in order to find {x'\e~P je \x). 
Here we explicitly evaluate Equation (6.29) for several values of M: 

M — 1: The expression for the thermal density matrix elements becomes 


{x'\p{p)\x) = 


{x'\e &\x) 

Tr( e -A^) 

(^%) 2ex P 


m(x) P n r 2 ■ ly' 2 \_ _m_( r r _ r x 2 l 

2 l2 A “ t_ 2 A ) 2 X) J 


(^Y fz> e *p\-^& 2+ 1* 2 )- 


i 


X\Ht V 2 


exp 


x + x'\ 2 ( fx-x') 2 4 + ( 2 


2x 


2x 


2( 


(6.30) 


where in the denominator we have set x = x' = x in order to evaluate the trace. 
We notice that this expression matches Equation (6.28) to first order in that 
is, Equation (6.30) is a high-temperature approximation of Equation (6.28). 
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M— 2: The expression for the thermal density matrix elements becomes 


(x'\p{P)\x) = 


(x'\e P*\x) 
Tr(e-P^) 



f 171 
U fr 2 P. 

In, exp 

[ [\x 2 + x 2 +\x’ 2 ) tfpKx 1 x) 2 + (x' Xi) 2 ] ] 

dxi 

"( 

m \ 
nh 2 p) 

/°oo eX P [ 

4 l3 [l x 2 + x 2 +lx 2 ) tfpliXi x) 2 + (x X!) 2 ]| 

cbqdx 


1 


la i 6 +< r 2 ) 

- 

' X + x' ’ 

2 (16 + ( 2 

' x - x' ’ 

2 8 + ( 2 

/ 4(8 + ( 2 ) c ^ p 

2x 

4 8 + ( 2 

2x 

4< 


(6.31) 


This is a slightly better approximation of Equation (6.28) for small (. 
M— 3: The expression for the thermal density matrix elements becomes 


(x'|pQ6)l*> 


(x'\e P*\x) 

Tr(tr^) 


3 

I 3 m Y r° nvn f mL ° 2 P 
{2nh 2 p) ^-oo eX P[ 6 

3 

( 3 m \ 2 roo l ma 2 p 
(27ih 2 l}) J-oo eX P[ 6 

1 / 243f(16 + ( 2 ) 

32(9 + C 2 )(27 + C 2 ) 


[\x 2 + x 2 + xf + \x' 2 ) - [(Xi -x) 2 + (x 2 -Xi ) 2 + (x'-x 2 ) 2 ]] dxidx 2 
(|x 2 + x 2 + x| + |x 2 ) - Jp^[(xi - x ) 2 + (x 2 - Xi ) 2 + (x- x 2 ) 2 ]| dxidx 2 dx 


- 

‘ x+ x!' 

2 f 27 + ( 2 

' x - x! ' 

2 (9 + C 2 )(36 + C 2 ) 

2x 

6 9-K 2 

2x 

6<T(27 + ( 2 ) 


(6.32) 


This is an even better approximation of Equation (6.28) for small (. 

M > 4: When you try to evaluate such integrals for a larger number M, in order to 
approach the limit M — > ■ oo, you will see that their evaluation takes a lot of 
computer power. They can be evaluated exactly in the present case of a har- 
monic oscillator; but in a more general case they cannot. 

We see from this series of explicit calculations that taking a finite value for M yields a 
high-temperature approximation of the thermal density matrix (approximately cor- 
rect for small (). The larger we choose M, the more the validity of the result extends 
to lower temperatures. 


6.3 Monte Carlo integration 

In Equation (6.13) and Equation (6.23) we have expressed matrix elements of the 
real- and imaginary-time propagators as integrals over many- dimensional spaces 
[for N particles and M time-slices, there are 3 N{M- 1) integration variables] . In this 
section we study a method for performing such integrals in practice, with a reason- 
able amount of computational power. 

6.3.1 one-dimensional uniform integration 

As a first example, we want to calculate a one-dimensional integral of the form 

J=f fix') dx, 

Jo 


(6.33) 
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where / : [0, 1] — *■ <D is an arbitrary function. 1 

Traditional Riemann integration of Equation (6.33) can, for example, be done 
with the rectangular rule: 



/= lim 

M—oo 


M -i 

y — x/ 


i 


rM f | 2_ 

^-m=l J M 


M 


= lim - 

M—oo 


M 


(6.34) 


If the function fix) is sufficiently well-behaved, the sum over the regular x-grid 
in Equation (6.34) can be replaced by a sum over a series of random numbers: if 
(xi,X 2 ,X 3 ,. ,.,xm) is a sequence of random numbers drawn independently and uni- 
formly from [0, 1], then (f(x m )) = fg f(x m )dx m and hence 


lim 

M—oo 


KLl /(*m) 

M 


{/(%)) = /. 


(6.35) 


This formulation is called a Monte-Carlo integral, after the city of Monte Carlo fa- 
mous for its gambling casinos. 

In Mathematica, we first define the function fix), for example 

In [378] : = f [x_] = x(l-x); 


and the correct answer for the integral, 

in [379] : = J = Integrate [f [x] , {x,0,l}] 
Out [379]= 1/6 


Since RandomReal [] generates random numbers drawn uniformly from [0,1], we 
calculate an average of M random numbers with 

in [380] : = Jmc [M_ Integer / j M> = 1] := Sum [f [RandomReal []], {M}] /M 


We can also simultaneously estimate the mean / = (fix)) and its standard error cr j ■ 

^1 {f 2 {x))-(f{xyf 


M 


from the same sequence (jci, Xz , . . . , xm)’ 


1 Choosing [0, 1] as the domain of integration is arbitrary; we could have chosen any real finite or infi- 

nite interval. 
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in[38i] : = Jmc [M_Integer/ ;M>=1] := Module [{x,fx}, 
(* the values x_m *) 
x = RandomReal [{0 , 1} , M] ; 

(* the values f(x_m) *) 
fx = f /8 x; 

(* return mean and standard error *) 
{Mean [f x] , Sqrt [Variance [f x] /M] }] 


For example, using M - 10000 we get a reasonable result for /: 


In [382] : = Jmc [10000] 

Dut[382] = {0. 166002, 0.000748734} 


6.3.2 one-dimensional integration with weight 


Next we wish to integrate a function / : [0, 1] — ► <D using a weight function p : [0, 1] — ► 
Rq satisfying fg p(x) dr = 1: 


7 = 


/' 


f(x)p(x) dx 


(6.36) 


In principle, we could define /(x) = f{x)p{x) and use the procedure of subsec- 
tion 6.3.1. In practice, there is a much more efficient procedure: we define the cu- 
mulative weight 

qix) = f p(y)dy, (6.37) 

Jo 

which satisfies q{ 0) = 0, q{ 1) = 1, and is monotonically increasing and therefore 
uniquely invertible on [0, 1], since q'[_x) — p{x) > 0. Hence, using the variable substi- 
tution z = q(x) we can re-express Equation (6.36) as 


7 = 


f f[q 1 {z)]dz = [ 
Jo Jo 


g(z)dz, 


(6.38) 


where we have defined g(z) = f[q 1 (z)]. Now we can use the procedure of subsec- 
tion 6.3.1 on this function g : [0, 1] — ► (D : 


/ = 


lim 

M— *■ oo 


Im=l StenJ 
M 


lim 

M—*oo 


M 


(6.39) 


where (zi, Z 2 , . . . , z,yj) is a sequence of random numbers drawn uniformly and inde- 
pendently from [0, 1]. We will show in an example that this is a much more efficient 
choice than using the f{x) defined above. 

Consider, for example, the weight function 

p{x) = 101 x x 100 (6.40) 


which is sharply concentrated around x = 1 but remains nonzero throughout (0,1]: 
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Defining the function f(x) = f(x) pix), as suggested above, means that when we use 
integration Equation (6.35) on / then more than 90% of the random numbers x m do 
not contribute significantly to the Monte Carlo estimate of /: 

In [383] : = p [x_] = 101 * X~100; 

in [384] : = J = Integrate [f [x] *p [x] , -fx,0,l}] 

Out [384] =10 1/10506 

in[385] : = Jmcl [M_Integer/ ;M>=1] := Module [{x , fpx} , 

(* the values x_m *) 
x = RandomReal [{0 , 1} , M] ; 

(* the values f (x_m) *p(x_m) *) 
fpx = f [#] p [#] & /@ x; 

(* return mean and standard error *) 

{Mean [fpx] , Sqrt [Variance [fpx] /M] }] 

In [386] : = Jmcl [10000] 

out [386] = {0.00948913, 0.000480376} 



We see that with 10 000 random numbers we got an estimate that has a relative pre- 
cision of about 5%. 

Now we calculate the cumulative weight 

qW= [ p(y)dy = x 101 , (6.41) 

Jo 

which in this case we can invert to q~ l (x) = x 1/101 . Using this, we calculate a second 
estimate of /: 

in[387] : = q[x_] = Integrate [p [y] , {y,0,x}] 

Out[387] = X~101 

in[388] : = Jmc2 [M_Integer/ ;M>=1] := Module [{z ,x,fx} , 

(* the values z_m *) 
z = RandomReal [{0 , 1} , M] ; 

(* the values x_m = Inverse [q] (z_m) *) 
x = z~ (1/101) ; 

(* the values f (x_m) = g(z_m) *) 
fx = f /<§ x; 

(* return mean and standard error *) 

{Mean [f x] , Sqrt [Variance [f x] /M] }] 
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In [389] : = Jmc2 [10000] 

out [389] = {0 . 00951386 , 0 . 0000907012} 


This time we get a relative precision of about 1% using the same number of random 
numbers. 

Given a sequence of random numbers (zi , zz , ■ ■ . , Zm ) drawn uniformly and inde- 
pendently from [0, 1 ] , we notice that the sequence (xi , X 2 , . . . , %) with x m - q \z 

m) 

is distributed according to p: 

in [390] :=z = RandomReal [{0 , 1}, 10000]; 
in[39i] :=Histogram[z, Automatic, "PDF"] 



In [392] : = Z = X~(l/101) ; 

in[393] :=Histogram [x , Automatic, "PDF"] 



Therefore, Equation (6.39) calculates an average of fix) using random values of x m 
drawn independently from the probability distribution p(x). 

6.3.3 the Metropolis-Hastings algorithm 

In our specific example p{x) = 101 x x 1 00 , drawing random numbers from this dis- 
tribution was relatively easy since the cumulative distribution q{x) - x 101 was an- 
alytically invertible. In more general cases, and in particular for multidimensional 
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probability distributions, this cannot be done and a different method for generating 
the random numbers x m must be sought. 

The simplest such algorithm is the Metropolis-Hastings algorithm, which gen- 
erates a sequence of correlated random numbers [x\,X 2 , ..., %) that are asymptot- 
ically (M —i • oo) distributed according to a probability density p{x). It works as fol- 
lows: 

1 . Start with a random starting value x \ , and set n- 1 . 

2. Propose a candidate for x n+ \ by drawing from a probability distribution n(x n >— * 

x n+ 1)- 

3. Calculate the acceptance ratio P(x n x n+ \) = min |l, ) ■ 

4. Choose a random number w n+ 1 from the uniform distribution over [0, 1): 

• If P{x n !-*• x n+ i) > w n + i, then we accept the move to x n+ \. 

• If P{x n !-*• x n+ i) < w n+ \, then we reject the move, and set x n+ \ - x n . 

5. Increment n and go back to step 2. 

Let’s do an example: as a candidate distribution we use 


' — ► x') 


2 ^ if|x-x'|<A, 
0 otherwise, 


(6.42) 


which is symmetric and therefore = 1 simplifies the acceptance ratio P{x n *— * 

x n +\)- 

in [394] :=next [x_ , d_] := Module [{y, P, w} , 

(* propose a new point *) 
y = x + RandomReal [{-d,d}] ; 

(* acceptance probability *) 

P = If [y<0 II y>l , 0, Min [1 , p[y]/p[x]]]; 

(* Metropolis-Hastings accept/reject *) 
w = RandomReal [] ; 

If[P > w, acc++; y, rej++; x] ] 

in[39B] :=MHchain [xl_?NumericQ , d_?NumericQ, M_Integer/ ;M>=1] : = 

Module [-Q , 

(* reset the acceptance/rejection counters *) 
acc = rej = 0; 

(* generate the chain of x values *) 

NestList [next [#,d]&, xl, M— 1] ] 


This procedure indeed generates a sequence {x\,X2 ,...,xm) distributed according 
to p(x) without making reference to the cumulative distribution q(x) or its inverse 
q~ l {x)\ 

in [396] :=X = MHchain[l, 0.015, 10000]; 
in[397] : = acc/(acc + rej) // N 
out [397] = 0.520152 


120 


CHAPTER 6. PATH-INTEGRAL METHODS 


Notice that we have picked a step size d = A = 0.015 such that the acceptance ratio 
is about 50%, i.e., about half of the proposed moves are accepted. 

in [398] :=P1 = Plot[p[x], {x, 0, 1}, PlotRange -> All]; 
in [399] :=P2 = Histogram [X , Automatic, "PDF"]; 

In [400] : = Show [P2 , P 1 ] 



X 


How does this work? Assume that the Metropolis-Hastings algorithm ultimately (in 
the limit M — ► oo) generates a set of values x m distributed according to a function 
s(x). This distribution function s(x) is therefore invariant under the Metropolis- 
Hastings algorithm, meaning that the detailed-balance condition s(x)tt{x x')P(x i— 
x') = s(x')ji(x' >->• x)P(x' >->• x) must be satisfied. Inserting the definition of P{x >-► x'), 


s(x)n(x^ x')min 


p(x')7r(x'~ x) \ 
p(x)7r(x'— x') ) 


s{x')ti{x' "->■ x) min 1, 


p{x)n(x >-► x') 
p{x')7i[x' >-<• x) 


(6.43) 


Since p and n are nonnegative, we can modify this to 


, J 

5(x)7T(X 1 — *■ X )- 


L [p(x)7r(x !-*• x'), p{x')n{x' >-<• x)] 
p(x)7l{x !-*• x') 

, , min [pfx'l^fx' >->■ x), w(x)7r(x >-► x')l 

= s(xV(x' — x) — (6.44) 


p{x')n(x' !-*• x) 

and, noticing that the minimum on both sides of this equation is the same, 


s(x) s(x') 
p{x) p(x') ' 


(6.45) 


The only way this equation can be satisfied for all (x, x') is if s(x) cx p(x). Since both 
s(x) and p(x) are normalized, we conclude that s(x) = p(x): the stationary distribu- 
tion of the Metropolis-Hastings algorithm is indeed p(x), as desired. 

What are such sequences (xi, X 2 , . . . , x ;V ;) good for? Since we know that the points 
in such a sequence are distributed according to p(x), we can now approximate inte- 
grals of the form of Equation (6.36) with 


f 


J= / f{x)p{x) dx- 


l M 

'■ i im 77 E 

M — oo M m=1 


(6.46) 
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6.4 Path-Integral Monte Carlo 


We can now combine the multi-dimensional integrals of Equation (6.13) and Equa- 
tion (6.20) with the stochastic integration method of section 6.3. 

We continue with the one-dimensional harmonic oscillator of subsection 6.2.1, 
in particular with Equation (6.29) for the matrix elements of the thermal density ma- 
trix. Comparing Equation (6.29) with Equation (6.36), we identify the weight func- 
tion 


p{x i,X 2 ,...,x M -\) °c exp 


mo) 2 /5 
2 M 


1 M—l i 

1 2 x -1 2 1 2 

-x 0 + £ X m + -X M 

^ m= 1 ^ j 


mM 
2 h 2 p 


M 

Y. i x m ~ x m- 1 ) 
m = 1 


(6.47) 


The goal of this section is to construct a sequence of paths whose elements are dis- 
tributed according to this weight function: as shown in the example of In [400] , 
the set shall contain more paths with high weight p{x \ , X 2 , . . . , %) and fewer paths 
with low weight. Notice that we need not be concerned with the pre-factor of p, as 
the Metropolis-Hastings algorithm will automatically find the correct normalization 
[see Equation (6.45)]. 

The Metropolis-Hastings algorithm can now be set up to work in the space of 
paths, that is, in the space of vectors x - (xj , x 2 , . . . , Xm- i), in the exact same way as 
we had set it up in subsection 6.3.3: 


1. Start with a random starting path x (1) , and set n = 1. A useful starting point 
would be the path that interpolates linearly between the fixed end points xo 
and xm, which is x$ — x o + [ml M){xm~ x o). 


2. Propose a candidate path x (n+1) by drawing from a probability distribution 


71 (X 


(n) 


;(«+ 1) 


). There are many ways of proposing new paths, and the effi- 


ciency of the stochastic integration will depend strongly on the choices made 
at this point. The simplest choice for finding a candidate is to select a ran- 
dom index pe ,M-\] and then add a random number Ax to x*"\ so that 


An+l) 


An) 


+ &m,[iAx. 


3. Calculate the acceptance ratio P{x [n) ~ x ( ' i+1) ) = min(l, )■ 

For the simple candidate mechanism above, the probability density is sym- 
metric, n{x' n ' > <-*■ x ! "~ l) ) = n(x" n+V) <-*■ x n> ), which simplifies the calculation of 
the acceptance ratio. 


4. Choose a random number w n+ i from the uniform distribution over [0, 1): 

• If P(x (,i) • 3c (,i+1) ) > w n+ 1 , then we accept the move to 3c ( ' i+1) . 

• If P(x^ x ln+1) ) < w n+ 1 , then we reject the move, and set - x*” 1 . 

5. Increment n and go back to step 2. 


We notice that, in the form presented here, the algorithm generates a sample of 
paths from Xq to x M that approximates the desired path weight function, Equa- 
tion (6.47), in the same way that in In [396] we had calculated a sample of values 
distributed according to the weight function given in Equation (6.40); but it does not 
yet give us an estimate for the density matrix element |xo). 

Let’s set up such a calculation in Mathematica. To simplify the notation, the ac- 
tion of a path x at inverse temperature p is expressed in terms of the dimensionless 


1 

2 

3 

4 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

1 

2 

3 

4 

5 

6 

7 

8 


122 


CHAPTER 6. PATH-INTEGRAL METHODS 


path coordinates x m = x m lx with the length scale x = 
inverse temperature ( = fiho) - 



and the dimensionless 


Six, 1 3) = 


ma> 2 p 
2 M 


M - 1 


n*0+ E 


X" + - X 


M 


m-l 

If 

M 


mM 9 

+ Wp^ 1 {Xm ~ Xm ~ l] 


1 M-l 
1 , V- ~2 


E 


1 


X" + — X' 


m=l 


M 


M M _ 2 

+ ~7T ^ (-^m — 

> m=l 


(6.48) 


With x-xlx- (jt 0 , xi,..., x M ) and z = ( we calculate this action: 


in[40i] reaction [x_/ ;VectorQ [x] &&Length [x] >=3 , z_] : = 

With [{M=Length [x] -1} , 

( (x [ [1] ] ~2/2+Sum [x [ [m] ] “2 , {m, 2 , M>] +x [ [M+l] ] ~2/2) * (z/M) 
+ Sum [ (x [ [m+l] ] -x [ [m] ] ) “2 , -Cm , 1 , M}] * (M/z) ) /2] 


Given a path x, we find the next path via the Metropolis-Hastings algorithm, using a 
random step of size d: 

in [ 402 ] :=next [x_/; VectorQ [x,NumericQ] &&Length[x] >=3 , 
z_?NumericQ, d_?NumericQ] := 

Module [{mu , dx , xn , S , Sn , P , w} , 

(* which point to modify *) 

(* (leave end points fixed!) *) 
mu = Randomlnteger [{2 .Length [x] -1}] ; 

(* by how much to move the point *) 
dx = RandomReal [{-d,d}] ; 

(* the new path *) 
xn = x ; xn [ [mu] ] += dx ; 

(* calculate path actions *) 

S = action[x,z] ; 

Sn = action[xn,z] ; 

(* acceptance probability *) 

P = Min [1 ,Exp [S-Sn] ] ; 

(* acceptance or rejection *) 
w = RandomReal [] ; 

If [P > w, acc++; xn, rej++; x] ] 


The Path-Integral Monte Carlo (PIMC) algorithm for generating a sample of u paths 
between i'o = xO and i.y; = xM taking M = M steps, at dimensionless inverse temper- 
ature ( = z, taking a random step Ax = d on average, looks thus: 

in [403] :=PIMCpaths [{xO_?NumericQ , xM_?NumericQ} , 

M_Integer/;M>=2, z_?NumericQ, d_?NumericQ, 
u_Integer/;u>=l] := 

Module [{x} , 

(* start with the straight path *) 
x = xO + Range [0,M]/M * (xM-xO) ; 

(* reset acceptance/rejection counters *) 
acc = rej = 0; 
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(* iterate the Metropolis-Hastings algorithm *) 
NestList [next [# , z , d] & , x, u]] 


Here is a graphical representation of 10 5 paths between in = 0 and % = 1 using M = 
20 imaginary- time slices, at a dimensionless inverse temperatures of £ = 10, 1,0.1,0.01: 


in [404] :=With[{z=l , M=20 , d=0 . 5 , u=10~5}, 
t = z * Range [0 , M] /M ; 
p = PIMCpaths [{0, 1}, M, z, d, u] ; 
DensityHistogram[Flatten [Transpose [{#,t}]&/@p, 1] , 
{Automatic, {-(z/(2M)), z(l+l/(2M)), z/M}}, 
{"Log", "PDF"}]] 


£=10 M=20 



X 


£=1 M=20 



X 


£=0.01 M=20 



-10 12 


x 


For smaller values of £, corresponding to higher temperatures, the paths are more 
and more concentrated around the straight path (the “least action” path of classical 
mechanics, indicated in red) . 

We notice that the output of PIMCpaths, which is a sequence of paths, does not 
directly allow us to calculate the matrix element (x'\e~^ je \x) from Equation (6.29); 
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instead, we can only evaluate integrals of the form of Equation (6.46). In what fol- 
lows, we will see what expectation values we can calculate directly from such path 
sequences. 


6.4.1 calculating the density 


The first observable quantity we wish to calculate is the thermal particle density 


p{x) = <x[p|x> = 


(. x\e P^\x) 
f (x'\e~^\x')dP N x' 


(6.49) 


Both the numerator and the denominator of this expression are diagonal matrix el- 
ements of e~ ^ and can be written as closed path integrals, where the end point is 
equal to the starting point of the paths: using Equation (6.20) with xq — xm - x and 

Xq — Xjyi — X , 


p{x) = lim 

M— *■ oo 


/exp 

[ ■ - 1 ( l m>) + LmZl V &m) + 5 V(.X M )) -§^Lm=l\\Xm-X m - 1 II 2 ] 

j3 N* 

d x ^ * ’ ’ d j 

/exp 

\-U\ v ^ + ^nl\ V& m )+ i Il3' m -2' m _ 1 ll 2 ] 

\d 3N x' 0 d 3N x[---d 3N x' M 


Notice that the denominator contains one more integration variable, d iN x' 0 , as re- 
quired in Equation (6.49). Assume now that we have an infinite sequence of closed 
paths (xm - x 0 ) through 3A/-dimensional configuration space, with asymptotic dis- 
tribution proportional to 


P{xq,xi,...,xm-i) oc exp 


P (l 1 

-V[x 0 )+ Y V(x m ) + -V{x M ) 

JVl Z yyj 1 Z 


mM 

2H 2 p 


M 


Y. II Xfn 

m= 1 


Xm - 1 


(6.51) 


2 


Of all these paths, the numerator of Equation (6.50) contains only those that start 
and end at x, while the denominator contains all of the paths: 


p(x) = 


number of closed paths starting from and ending in x 
number of closed paths 


(6.52) 


We notice further that since these paths are closed, we cannot tell which point is 
their starting point and which is their end point; this insight improves the statistics 
of Equation (6.52) by a factor of M to 


p(x) = 


number of closed paths containing x 
M times the number of closed paths 


(6.53) 


In practice, calculating the density in 3 ]V-dimensional configuration space thus boils 
down to making a 3/V-dimensional histogram of all the points contained in all the 
closed paths of the sequence. We will illustrate this with our harmonic oscillator 
example. 


thermal density of a harmonic oscillator 


The thermal density of a harmonic oscillator can be calculated analytically from 
Equation (6.28): 


p{x) = (x\ p{p)\x) = 


exp[-(|) 2 tanh(|)] 


(6.54) 
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We calculate a sequence of closed paths, where the beginning/end of the path is mo- 
bile as well, through a slight modification of In [402] and In [403] : the last element 
of the list x = (xo, X| , . . . , xm- t) has been chopped off, since it is identical to the first 
element. However, if we re-use the code of In [402] with this modification, we no- 
tice quickly that the convergence of the path sequence to the desired distribution is 
terribly slow. The reason is apparent: if we move only one point of the path at a time, 
it takes a long time until the entire path can move to a different place. The scale of 
motion of a single point on the path is given by the thermal de Broglie wavelength 
of the particle, and therefore goes to zero at high temperature; at the same time, the 
scale of motion of the entire path is given by the thermal width of the density, which 
becomes larger at high temperature. We must therefore introduce a second type of 
move, one that displaces the entire ring. 

The first type of move remains the same: displace one point on the path by a 
random distance. 

in [405] :=nextCl [x_/ ; VectorQ [x ,NumericQ] &&Length [x] >=2 , 
z_?NumericQ, d_?NumericQ] := 

Module [{mu , dx , xn , S , Sn , P , w)- , 

(* which point to modify *) 

mu = Randomlnteger [{1 , Length[x] }] ; 

(* by how much to move the point *) 
dx = RandomReal [{-d, d}] ; 

(* the new path *) 
xn = x ; xn [ [mu] ] += dx ; 

(* calculate path actions *) 

S = act ion [Append [x, First [x]], z] ; 

Sn = action [Append [xn, First[xn]], z] ; 

(* acceptance probability *) 

P = Min[l, Exp[S-Sn]]; 

(* acceptance or rejection *) 
w = RandomReal [] ; 

If [P>w, accl++;xn, rejl++;x]] 


The second type of move displaces the entire path (ring) by a random distance: 

in [406] :=nextC2 [x_/ ; VectorQ [x ,NumericQ] &&Length [x] >=2 , 
z_?NumericQ, d_?NumericQ] := 

Module [{dx , xn , S , Sn , P , w} , 

(* by how much to move the points *) 
dx = RandomReal [{-d, d}] ; 

(* the new path *) 
xn = x + dx ; 

(* calculate path actions *) 

S = act ion [Append [x, First [x]], z] ; 

Sn = action [Append [xn, First [xn]], z] ; 

(* acceptance probability *) 

P = Min[l, Exp[S-Sn]]; 

(* acceptance or rejection *) 
w = RandomReal [] ; 

If [P>w, acc2++;xn, rej2++;x]] 
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At every iteration, we choose a move of type 1 with probability f and a move of type 
2 with probability 1 - f : 

in [ 407 ] :=nextC [x_/ ; VectorQ [x ,NumericQ] &&Length [x] >=2 , 

z_?NumericQ, dl_?NumericQ , d2_?NumericQ , 
f _?NumericQ] := 

If [RandomReal [] >f , nextC2[x,z,d2] , nextCl [x,z,dl]] 


Construct a sequence of closed paths: 

in [408] :=PIMCpathsC [xO_?NumericQ , 

M_Integer/;M>=2, z_?NumericQ, dl_?NumericQ , 
d2_?NumericQ , f_?NumericQ, u_Integer/ ;u>=l] := 

Module [{x> , 

(* start with the trivial path at xO *) 
x = Table [xO, {M}] ; 

(* reset acceptance/rejection counters *) 
accl = rejl = acc2 = rej2 = 0; 

(* iterate the Metropolis-Hastings algorithm *) 
NestList [nextC [# , z, dl ,d2 , f ]& , x, u] ] 


The density is found by plotting a histogram of all the points contained in all the 
paths: 

in[409]:=With[{z=l, M=20 , dl=0.45, d2=3, f=0.5, u=10~5}, 

X = PIMCpathsC [0 , M, z, dl, d2 , f, u] ; 

Histogram [Flatten [X] , Automatic, "PDF"] 


£=10 M=20 


£= 1 M=20 


X 


X 



We see that these densities match the analytic expressions [red lines, Equation (6.54)] 
within statistical uncertainties. While the spread of each path decreases with de- 
creasing £ (see In [404] ), the overall size of the density profile increases with de- 
creasing £; hence the need for two different kinds of random moves above. We can 
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see this decreasing ring size by moving each ring such that its center of gravity is at 
x = 0, and plotting a histogram of the resulting points: 

i in[4io] :=Histogram[Flatten [#-Mean [#] &/@X] , Automatic, "PDF"] 
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In the classical limit (£ — «■ 0), the closed paths are reduced to loops of zero length, i.e. t 
points, but these points are distributed over a large region in space. This is a simple 
example of how path integrals provide an intuitive picture of the classical limit of 
quantum mechanics. 
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eigenvalues, 27, 48, 50 
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